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disease  cap  be  identified,  the  smaller  subset  (perhaps  200-500  per  10,000)  at  high  risk 
of  sudden  d^ath  can  then  be  removed  from  flying  status.  The  potential  savings  would  clearly 
justify  a  systematic  exploration  of  the  possibility  that  a  cost-effective  system  for  routine 
electrocardiographic  (ECG)  mapping  at  rest  and  with  exercise  can  be  developed  that  would 
identify  95%  or  more  of  these  pilots.  Recent  findings  from  a  number  of  centers  and  analysis 
of  available  data  suggest  that  an  overall  sensitivity  for  detecting  ischemic  coronary  diseast 
of  95%-98%  may  be  possible  along  with  a  specificity  in  the  same  range. 

This  technical  report  details  previous  work  by  us  and  a  number  of  investigators 
defining  the  anatomic  and  physiologic  basis  for  the  simulation  of  depolarization/ 
repolarization  of  the  heart,  i.e.,  QRS-T  changes  on  the  electrocardiogram  (ECG).  A  detailed 
design  of  new  programs  to  expand  our  fine  grid  (1  mm3)  propagation  simulation  of  depolar¬ 
ization  (QRS)  to  include  repolarization  (ST-T)  is  presented.  This  forward  computer  model 
will  include  a  family  of  inhomogeneous  torsos  representative  of  Air  Force  pilots  at  rest 
and  at  exercise.  Using  Bayesian  Statistics  and  this  model,  we  expect  to  develop  criteria 
for  local,  body  surface  ECG  changes  that  are  specific  for  each  of  2  (or  4)  main  coronary 
arteries  and  quantitative  for  the  degree  of  ischemia  or  infarction  or  both  in  each  of  these 
coronary  distributions.  An  optimal  subset  (8-32  leads)  of  a  176-lead  rest  and  exercise 
ECG  map  will  be  chosen  that  gives  the  greatest  sensitivity,  specificity,  and  diagnostic 
accuracy  in  the  detection  of  myocardial  infarcts  and  significant  ischemia  from  coronary 
artery  disease.  When  validated,  this  database  will  be  used  to  develop  a  portable,  multi¬ 
lead  ECG  Mapping  System.  This  high-fidelity,  automated,  digital  system  can  be  used  in 
flight  surgeon's  offices  worldwide  for  resting  ECGs  and  exercise  testing  to  improve  sig¬ 
nificantly  Lhe  detection  of  regional  coronary  iscemia  and  quantitate  infarct  size  in  the 
asymptomatic  pilot.  We  expect  to  achieve  a  predictive  accuracy  of  95%  or  greater. 
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EXECUTIVE  SUMMARY 


Ten  to  20%  of  all  apparently  healthy  people  over  35  have  significant  but 
undiagnosed  coronary  artery  disease.  For  15  to  30%  of  these  people,  sudden 
death  is  the  first  symptom  of  heart  disease  and  as  many  as  half  of  those  who 
die  have  had  a  prior  unrecognized  infarct  (heart  attack).  Applying  these 
statistics  to  pilots,  at  least  1  pilot  per  10,000  over  the  age  of  35  who  fly 
an  annual  average  of  300  hours  in  high-performance  aircraft  will  die  each 
year  at  the  controls  resulting  in  loss  of  the  aircraft.  Based  on  published 
reports,  the  Air  Force  loses  from  20  million  to  as  much  as  50  million  dollars 
per  year  from  aircraft  lost  this  way.  Furthermore,  it  follows  that  there  are 
10-15  unrecognized  new  infarcts  per  year  per  10,000  of  these  same  flyers. 
Improving  the  identification  of  coronary  disease  in  these  pilots  would  lower 
the  number  flying  with  unrecognized  heart  disease.  This  would  decrease  the 
chances  of  sudden  death  of  pilots  at  the  controls  and  reduce  the  number  of 
aircraft  lost. 

This  report  describes  proposed  enhancement  of  an  existing  computer  model 
to  improve  ECG  criteria  and  to  define  an  optimal  set  of  ECG  electrodes  for 
the  identification  of  significant  coronary  artery  disease  with  and  without 
myocardial  infarction.  The  ultimate  goal  is  the  development  of  a  portable 
ECG  Mapping  Cart  that  can  be  used  in  flight  surgeon's  offices  worldwide  for 
resting  and  exercise  ECGs  to  improve  coronary  artery  disease  detection  and 
quantitation  of  infarct  size  in  the  asymptomatic  pilot. 
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OPTIMAL  ECG  ELECTRODE  SITES  AND  CRITERIA  FOR  DETECTION 
OF  ASYMPTOMATIC  CORONARY  ARTERY  DISEASE 
AT  REST  AND  WITH  EXERCISE 


OVERVIEW  AND  OBJECTIVES 

(University  of  Southern  California  Perspective) 


Role  of  Model  Building 

Historically,  scientific  models  have  been  used  to  explain  inherently 
complex  natural  systems.  In  biology,  specifically  in  electrocardiography  and 
electrophysiology,  Einthoven,  the  father  of  electrocardiography,  viewed  the 
heart  as  a  single  fixed  locus  dipole.  This  simple  model  of  a  complex  system 
of  cellular  electrical  generators  still  explains  many  aspects  of  the  heart  as 
an  electric  current  source.  Einthoven  and  associates  used  an  equilateral 
triangle  as  a  model  of  the  human  torso  to  relate  the  electrocardiographic 
signals  from  the  three  standard  limb  leads  (1).  This  model  explains  so  many 
aspects  of  the  complex,  inhomogeneous,  irregularly  shaped  human  torso  that  it 
forms  a  common  basis  for  all  clinical  electrocardiographers  to  define  speci¬ 
fic  changes  seen  on  ECGs.  The  concept  of  a  mean  electrical  axis  for  the  ECG 
depends  on  both  the  equivalent  dipole  model  of  the  heart  as  a  current  source 
and  the  equilateral  triangle,  or  hexagonal  reference  system,  as  a  model  of 
the  complex,  inhomogeneous  volume  conductor  through  which  these  currents 
travel . 


Models  of  Ventricular  Depolarization 

Sir  Thomas  Lewis,  across  the  English  Channel  from  Einthoven,  went  to  the 
source  and  made  elegant  measurements  on  the  heart  itself  using  Einthoven's 
string  galvanometer  (2-4).  From  these  measurements,  he  developed  the  dipole 
layer  model  to  explain  the  propagating  wave  of  electrical  excitation  he  found 
on  the  epicardium.  For  years,  his  model  of  intramural  excitation  was  used  to 
explain  the  complex  waveforms  of  conventional  ECGs  and  of  the  ECGs  of  pa¬ 
tients  with  myocardial  infarction,  bundle  branch  blocks,  and  hypertrophy  of 
either  ventricle.  This  early  model  of  excitation  of  the  ventricles  was  re¬ 
vised  considerably,  particularly  the  activation  of  the  septum,  after  Allan 
Scher  made  his  classic  studies  on  intramural  excitation  using  small,  intra¬ 
mural  plunge  electrodes  (5-8).  Durrer  and  associates  reported  similar  find¬ 
ings  (9-13).  However,  the  uniform  propagating  dipole  layer  was  retained  as 
the  model  of  the  excitation  wave  fronts  moving  through  the  heart. 

When  we  developed  the  computer  simulation  of  excitation  of  the  human  ven¬ 
tricle,  we  decided  to  build  the  simulation  using  measured  anatomy,  resistiv¬ 
ities,  and  electrophysiology.  From  the  intramural  electrograms  we  sought 
data  to  define  the  dipole  layer  separation  and  the  source  strength  of  the 
positive  and  negative  poles  of  the  uniform  dipole  layer.  We  discovered  that 
the  dipole  layer  model  of  the  active  electromotive  surface  was  inconsistent 
with  measured  data  (14),  and  a  new  model  of  the  propagating  wave  front  needed 


to  be  developed  from  the  measured  intramural  data.  Using  a  fine  grid  (1  mm3) 
digital  computer  simulation  of  the  propagating  wave  of  depolarization  (QRS 
changes)  of  the  human  heart  we  found  that  small  myocardial  infarcts  no  larger 
than  1-1.5  cm  in  diameter  and  0.4-0. 5  cm  thick  routinely  produced  detectable 
changes  in  the  simulated  ECGs  (15,16).  In  clinical  studies,  we  confirmed 
that  small,  clinically  "silent"  infarcts  can  be  routinely  detected  from  QRS 
changes  in  good  quality  resting  electrocardiograms  and  vectorcardiograms 
(VCGs)  when  prior  tracings  are  available  for  comparison.  We  found  that 
serial  change,  or  lack  of  it,  in  high-gain,  high-fidelity  ECG/VCG  is  predic¬ 
tive  of  serial  ventriculographic  change,  or  lack  of  it,  with  98%  accuracy 
(17). 

We  expect  this  to  have  particular  relevance  to  the  pilot  population  in 
whom,  for  reasons  discussed  in  the  Air  Force  perspective,  it  can  be  assumed 
that  more  than  half  of  the  myocardial  infarcts  suffered  will  have  been  clini¬ 
cally  unrecognized,  or  "silent".  Our  model  incorporated  measured  anatomy  and 
physiology  of  ventricular  excitation  that  allowed  us  to  develop  and  validate 
new  QRS  criteria  for  quantitating  myocardial  infarct  size  (18-23).  This 
success  has  given  us  confidence  that  the  same  strategy  extended  to  include 
repolarization  in  the  fine  grid  model  is  likely  to  produce  more  sensitive  and 
specific  ST-T  criteria  for  detecting  ischemia  due  to  coronary  disease. 


Models  of  Ventricular  Repolarization 

Most,  if  not  all,  models  of  ventricular  repolarization  are  indebted  to 
some  extent  to  the  ventricular  gradient  model  of  Wilson  (24,25)  and  the  many 
studies  refining  it  (26-49).  Harumi,  Burgess,  and  Abildskov  were  the  first 
to  develop  a  model  of  the  T  wave  (50).  Thiry  and  Rosenberg  later  developed  a 
simulation  of  repolariza+ion  in  3  dimensions  (51).  Miller  and  Geselowitz 
developed  a  4  mm3  grid  heart  model  in  a  homogeneous  torso  in  which  they 
manually  entered  changes  in  excitation  and  recovery  parameters  (52-53). 
Essentially,  we  will  be  automating  the  grid  heart  model  of  Miller  and 
Geselowitz  in  a  1  mm3  grid  heart  in  an  inhomogeneous  torso.  The  success  of 
these  models  in  simulating  typical  12-lead  ECG  changes  of  ischemia,  injury, 
and  infarction  by  incorporating  action  potential  waveshapes  and  distributions 
typical  of  those  measured  in  the  laboratory  makes  the  likelihood  of  success 
with  the  fine  grid  version  quite  high. 


The  Role  of  the  Computer  Simulation  of  QRS-T 
in  the  Development  of  Quantitative  Criteria 
and  Optimal  Lead  Sets 


Simulation  of  Body  Surface  Maps  at  Rest 


Briefly,  the  effort  is  to  develop  a  mathematical  generator  for  repolari¬ 
zation  and  simulate  the  myocardial  depolarization/recovery  process.  Our 
method  of  constructing  the  repolarization  generator  makes  use  of  an  equation 
containing  three  arbitrary  parameters  that  assume  values  which  permit  the 
time  course  of  the  EMF  (electromotive  force  or  electromotive  surface  n  the 
heart)  across  the  cell  membrane  to  be  approximated.  The  three  parameters 
characterizing  the  action  potential  are  the  recovery  potential,  rate  of 


repolarization,  and  the  functional  refractory  period.  With  these  3  para¬ 
meters  and  4  constants  (initial  resting  potential,  time  to  depolarize,  peak 
depolarization  potential,  and  a  time  constant  which  represents  the  time  for 
the  local  region  to  reach  equilibrium),  each  discrete  point  within  the 
myocardium  can  have  a  unique  waveform  that  represents  it. 

This  equivalent  generator  is  then  combined  with  a  Gelernter-Swihart  (GS) 
solution  for  the  heart-to-chest  transfer  impedances.  Combining  the  generator 
and  transfer  impedances  yields  a  mathematical -physical  simulation  containing 
arbitrary  parameters  that  can  be  adjusted  such  that  numerically  computed  body 
surface  T-wave  distributions  can  be  made  to  agree  with  those  measured  on 
human  subjects.  Once  a  verified  model  is  developed  and  fitted  to  resting  and 
exercise  normal  map  data,  then  a  minimum  set  of  chest  electrodes  can  be 
selected  which  optimizes  specificity  and  sensitivity  of  the  model  to  ischemia 
with  respect  to  ECG  waveform  changes. 


Simulation  of  Exercise  Maps 

It  is  anticipated  that  the  coefficients  and  constants  in  the  model  for 
simulating  ECGs  will  have  to  be  altered  to  account  for  changes  that  occur 
with  exercise.  Subjects  during  stress  testing  usually  exercise  in  an  upright 
position;  their  torso  geometry  and  resistivities  change  as  they  increase 
their  tidal  volumes  and  blood  flow  above  resting;  their  heart  position 
changes  as  does  the  distance  between  their  surface  ECG  leads  and  local  car¬ 
diac  generators.  It  follows  that  transfer  impedances  from  local  cardiac 
generators  to  body  surface  leads  may  be  significantly  changed  between  quiet 
resting  and  exercise  since  they  vary  inversely  as  the  square  of  the  distance 
between  generator  and  surface  electrode.  Trans-thoracic  resistivity  measure¬ 
ments,  measured  thoracic  geometry,  and  heart  location  by  2-D  echocardiograms 
during  quiet  supine  rest  and  various  levels  of  upright  exercise  will  be  used 
to  adjust  the  torso  simulation  with  the  Gelernter-Swihart  program  to  compute 
a  series  of  transfer  impedances  between  heart  and  surface  ECG  leads  for  these 
various  intensities  of  exercise. 

Additional  changes  in  the  repolarization  generator  parameters,  such  as 
recovery  rate  and  functional  refractory  period,  will  be  introduced.  In¬ 
creased  heart  rate  and  heart  muscle  temperature  are  known  to  change  the  shape 
of  the  monophasic  action  potential  during  recovery.  A  set  of  normal  human 
subjects  will  be  employed  to  adjust  the  G-S  transfer  impedances  at  exercise 
to  bring  simulated  QRS-T  wave  potential  distributions  at  exercise  into 
agreement  with  observed  data. 

Our  methods  of  cycle  selecting  and  averaging  exercise  body  surface  ECGs 
will  be  the  same  as  at  rest  since  the  principles  are  invariant  with  respect 
to  differences  between  rest  and  exercise.  Adjusting  the  simulation  to  ac¬ 
count  for  exercise  QRS-T  wave  potential  distributions  will  be  accomplished  by 
the  method  of  residuals.  Our  aim  will  be  to  reduce  the  residuals  at  exercise 
to  a  level  commensurate  with  those  attained  from  resting  QRS-T  wave  distribu¬ 
tion.  This  aim  comes  out  of  the  desire  to  maintain  the  same  level  of  cer¬ 
tainty  in  the  simulation  of  rest  and  exercise  ECGs.  When  similar  confidence 
limits  exist  on  both  sets  of  data,  this  simplifies  the  interpretation  of  the 
repolarization  generator  differences  between  rest  and  exercise. 
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Simulation  Model  of  QRS-T  in  Preliminary  Hypothesis  Testin 

A  computerized  body  surface  ECG  map  simulation  will  provide  an  alterna¬ 
tive  method  of  studying  pathophysiological  mechanisms  and  generating  less 
commonly  observed  body  surface  maps.  A  matrix  transformation  between  differ¬ 
ence  properties  (the  difference  between  multilead  ECG  maps  at  rest  and  exer¬ 
cise)  and  the  severity  of  coronary  occlusion  will  form  the  basis  for  the 
analysis  of  regional  ischemia  in  a  selected  patient  population.  This  matrix 
transformation  takes  the  form: 

4  E  =  T-C  +  £ 

where  the  elements  of  the  vector  4E  are  the  changes  in  ECG  voltages  between 
rest  and  exercise.  The  vector  C  is  composed  of  elements  whose  values  are  the 
regional  degree  of  coronary  occlusion  from  patients  with  single  vessel  di¬ 
sease  documented  by  coronary  angiograms.  The  unknown  T  matrix  of  weighting 
factors  will  be  formed  from  this  data.  The  validated  ECG  simulation  will 
provide  data  for  C  andAE  vectors  not  available  in  the  learning  set  and  will 
be  used  to  fill  out  the  unavailable  elements  in  the  T  array  of  weighting 
factors.  The  specific  hypothesis  of  this  simulation  effort  is  that  exercise- 
induced  ischemia  and  QRS-T  changes  in  body  surface  maps  can  be  transformed 
into  degrees  of  regional  coronary  artery  disease.  This  hypothesis  is  cast  as 
mathematical  transformations  between  exercise-induced  QRS-T  changes  in  spe¬ 
cific  regions  of  the  torso  and  the  ischemia  produced  in  regional  myocardium 
by  specific  coronary  disease. 


OVERVIEW  AND  OBJECTIVES 
(U.S.  Air  Force  Perspective) 
by 

Gil  D.  Tolan,  Colonel,  USAF,  MC 

(Author's  note:  this  summary  from  an  earlier  statement  of  work  is  presented 
here  because  it  gives  those  from  outside  the  Air  Force,  including  our 
consultants,  a  perspective  from  "inside"  that  would  not  otherwise  be 
available. ) 

Although  coronary  artery  disease  mortality  in  the  United  States  decreased 
22.6%  between  1969  and  1977,  at  670  deaths  per  100,000  men  per  year,  coronary 
artery  disease  still  remains  the  number  one  cause  of  death  in  this  country. 
Coronary  artery  disease  is  the  cause  of  70%  of  the  nontraumatic  deaths  in  the 
active  duty  U.S.  Air  Force.  Rated  flying  personnel  are  not  immune  to  devel¬ 
oping  coronary  artery  disease.  Each  year  approximately  50  rated  pilots  and 
navigators  on  flying  status  are  admitted  to  USAF  hospitals  where  the  dis¬ 
charge  diagnosis  is  coronary  artery  disease  (ICDA  codes  410  through  414). 
Each  year  approximately  8  rated  personnel  on  flying  status  suffer  fatal 
myocardial  infarctions.  It  is  hard  to  say  how  many  USAF  aircraft  accidents 
per  year  are  due  to  coronary  artery  disease  because  frequently  there  is  no 
heart  found  to  autopsy  after  the  crash,  and  because  the  on-board  recorder 
does  not  record  the  ECG.  Even  if  a  pilot  suffered  an  acute  myocardial 
infarction  in  flight,  there  would  not  be  sufficient  time  between  that  event 
and  when  the  rest  of  the  heart  died  in  the  crash  to  distinguish  the  acute 
infarction.  It  takes  hours  to  days  for  the  necrosis  of  an  infarction  to 
become  evident  at  autopsy.  The  Armed  Forces  Institute  of  Pathology  has  a 
policy  of  not  attributing  an  aircraft  crash  to  coronary  artery  disease  merely 
because  an  old  myocardial  infarction  scar  is  found  or  because  coronary  ather¬ 
osclerosis  is  present.  Today  44%  of  the  rated  force  is  over  35,  the  age 
susceptible  for  coronary  artery  disease.  Coronary  artery  disease  is  not 
confined  to  just  older  pilots.  Autopsy  studies  on  military  men  killed  in 
combat  in  both  Viet  Nam  and  Korea  found  coronary  artery  disease  in  many  young 
troops.  Although  we  do  not  know  how  many  aircraft  accidents  per  year  are  due 
to  coronary  artery  disease,  we  do  know  that  approximately  six  times  per  year 
a  USAF  aircraft  for  no  apparent  reason  impacts  the  ground  with  no  prior 
communication  from  the  pilot  and  no  attempt  to  eject.  Some  suspect  that  the 
pilot  was  trying  to  recover  an  out-of-control  aircraft  right  up  to  the  end, 
but  others  believe  the  lack  of  any  communication  of  trouble  implies  sudden 
human  incapacitation.  Although  there  are  many  causes  of  sudden  incapacita¬ 
tion,  common  things  occur  commonly,  and  coronary  artery  disease  is  the  most 
common  cause  of  sudden  incapacitation  in  American  males. 

Today  the  screening  for  coronary  artery  disease  in  the  USAF  rated  force 
relies  solely  on  the  standard  12-lead  ECG.  The  Framingham  14-year  follow-up 
of  5000  men  and  women  considered  free  of  heart  disease  at  onset  found  that 
the  standard  12-lead  ECG  was  normal  in  2  out  of  every  3  cases  which  suffered 
a  myocardial  infarction  in  the  subsequent  2  years.  This  number  represents  a 
33%  sensitivity  for  the  standard  12-lead  ECG  as  a  predictor  of  a  myocardial 
infarction  in  the  next  2  years.  The  Framingham  study  also  found  that  in 
those  cases  with  scars  from  old  myocardial  infarctions  found  at  autopsy,  36% 


of  these  were  unrecognized  during  life  by  both  the  patient  and  his  doctor. 
In  16%  of  all  coronary  cases,  the  first  symptom  was  sudden  death. 

Most  USAF  flyers  enjoy  flying;  they  get  paid  extra  to  be  on  flying 
status,  and  often  their  career  progression  requires  that  they  be  fit  to  fly. 
Although  most  flyers  have  IQ's  three  standard  deviations  above  the  mean,  they 
nearly  all  score  three  standard  deviations  above  the  mean  on  the  MMPI  denial 
scale.  Their  ability  to  deny  risk  is  a  life-preserving  attribute  that  pre¬ 
vents  them  from  panicking  when  exposed  to  enemy  fire.  These  same  factors, 
however,  decrease  a  flyer's  propensity  to  report  his  chest  pain  to  his  flight 
surgeon. 

In  the  new  generation  high-performance  aircraft,  flyers  are  frequently 
exposed  to  greater  than  7  +G  forces.  By  the  end  of  this  decade,  two-thirds 
of  USAF  active-duty  pilots  will  be  flying  high-performance  aircraft.  These 
aircraft  place  unprecedented  demands  on  the  heart  to  keep  the  brain  perfused. 
Within  30  seconds  of  pulling  7  +G's,  the  pooling  of  blood  in  the  lower  lobes 
of  the  lungs  causes  the  alveolae  to  collapse,  creating  a  ventilation  perfu¬ 
sion  abnormality  such  that  in  healthy  nonsmokers  their  mean  arterial  oxygen 
tension  drops  from  91  to  46  mmHg.  At  the  very  time  the  demand  for  cardiac 
output  is  maximum,  the  oxygen  supply  to  the  heart  is  being  decreased  by  both 
this  ventilation  perfusion  abnormality  and  by  the  shortening  of  diastolic 
filling  time  from  the  faster  heart  rate.  Since  flow  of  a  liquid  through  a 
tube  is  inversely  proportional  to  the  fourth  power  of  the  radius,  a  slight 
decrease  in  a  coronary  artery's  diameter  can  have  a  big  effect  on  the  amount 
of  blood  that  gets  by  that  obstruction  should  the  flow  remain  laminar. 
Because  a  slight  obstruction  can  cause  turbulence,  a  slight  coronary  obstruc¬ 
tion  can  be  more  hemodynamical ly  significant  than  implied  by  the  percent 
diameter  narrowing.  Furthermore,  a  healthy  coronary  artery  can  increase  flow 
ninefold  by  dilating  during  peak  workloads.  By  preventing  dilation,  circum¬ 
ferential  atherosclerosis,  even  if  minimal  on  angiography,  can  profoundly 
reduce  myocardial  oxygen  supply  during  peak  loads.  Even  without  causing  a 
myocardial  infarction,  an  asymptomatic  coronary  narrowing  can  decrease  maxi¬ 
mum  cardiac  output  during  sustained  high  +G  maneuvers.  If  during  a  high  +G 
maneuver  the  pilot's  brain  is  not  perfused,  then  under  the  best  of  circum¬ 
stances,  it  will  take  him  on  the  average  15  seconds  to  recover  from  his  loss 
of  consciousness.  By  then  it  may  be  too  late  to  recover  his  supersonic 
aircraft.  Thus,  a  slight  coronary  artery  narrowing  which  might  never  cause 
an  infarction  on  the  ground  could  prove  fatal  under  the  sustained  high  +G 
profiles  of  the  new-generation  aircraft. 

For  all  the  reasons  cited  above,  flying  safety  requires  that  we  not  wait 
for  the  flyers  to  come  to  their  flight  surgeon  with  complaints  of  chest  pain 
before  we  begin  worrying  whether  they  have  coronary  artery  disease.  We  need 
to  be  able  to  detect  coronary  artery  disease  while  it  is  still  asymptomatic. 
At  this  time  the  National  Institutes  of  Health  have  no  plans  to  fund  research 
into  detecting  asymptomatic  coronary  artery  disease  because  "the  only  treat¬ 
ment  for  it  is  advice  to  modify  risk  factors  and  that  advice  should  be  given 
regardless  of  whether  coronary  artery  obstructions  already  exist."  If  medi¬ 
cal  technology  is  to  keep  pace  with  the  demands  of  new  aircraft  technology, 
then  it  is  up  to  the  U.S.  Air  Force  to  develop  a  better  system  to  detect 
coronary  artery  disease. 
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The  goal  of  this  contract  is  to  develop  a  system  which  is  better  than 
the  existing  method  at  distinguishing  rated  flyers  with  aeromedical ly  dis¬ 
qualifying  coronary  artery  disease  from  healthy  flyers.  This  system  includes 
detecting  both  coronary  artery  obstructions  (ischemia)  and  myocardial  infarc¬ 
tions  (scars).  Correctly  localizing  infarctions  and  coronary  obstructions 
with  a  screening  test  in  the  field  is  not  a  requirement  of  this  contract. 

Good  skin-to-electrode  contact  is  essential  to  accurate  electrocardio¬ 
graphy,  but  good  contact  during  stress  testing  is  difficult  to  maintain 
because  the  adhesive  may  loosen  as  the  patient  moves  and  perspires.  The  more 
electrodes  one  adds,  the  greater  this  problem.  Furthermore,  additional  elec¬ 
trodes  require  additional  time  to  apply,  more  channels  for  recording  and  more 
complex  analysis,  all  of  which  increase  the  cost.  For  these  reasons,  the 
U.S.  Air  Force  needs  to  use  enough  electrodes  in  both  resting  and  stress 
electrocardiography  to  ensure  flying  safety  without  adding  any  extra  elec¬ 
trodes  merely  to  determine  some  day  if  they  might  be  useful.  Any  addition  of 
electrodes  will  incur  a  considerable  expense  not  only  to  acquire  new  ECG 
carts  for  the  field,  but  also  to  train  USAF  physicians  to  interpret  the  new 
leads.  For  this  reason,  any  new  lead  must  be  based  upon  strong  scientific 
evidence  that  it  adds  important  information  not  currently  available.  Because 
the  expense  and  risks  of  both  false  positive  and  false  negative  errors  are 
considerable,  the  process  in  developing  such  scientific  evidence  must  be 
meticulous  in  every  detail.  Because  current  policy  requires  that  rated 
flyers  with  abnormal  stress  tests  be  grounded  until  proven  healthy  by  cardiac 
catheterization,  it  is  essential  that  any  new  screening  system  maintain  a 
specificity  of  95%  or  better. 

In  1974,  Victor  Froelicher,  Jr.,  Ma j ,  USAF,  MC,  reported  in  the  American 
Journal  of  Cardiology  on  a  group  of  1,390  asymptomatic  men  screened  at  the 
USAF  School  of  Aerospace  Medicine  for  latent  coronary  artery  disease  with  a 
double  Master's  two-step  test  followed  by  a  symptom-limited  treadmill  exer¬ 
cise  tolerance  test.  The  follow-up  period  ranged  from  4.1  to  8.4  years  (mean 
6.3).  End  points  for  coronary  artery  disease  were  angina  pectoris,  acute 
myocardial  infarction,  and  sudden  death.  Only  lead  X  (CC5)  was  always  recor¬ 
ded  during  exercise.  Leads  II,  III,  aVF,  V2,  and  V5  were  recorded  before 
exercise  and  during  recovery.  In  this  follow-up  study,  the  Master's  test  had 
a  sensitivity  of  61%  and  a  specificity  of  92%.  Subsequently,  this  institu¬ 
tion  has  recorded  four  leads  (X,  Yh,  Z,  and  CM5)  during  all  phases  of  the 
treadmill  test.  No  follow-up  study  of  this  lead  system  has  been  done,  but  we 
know  that  60%  of  the  asymptomatic  flyers  with  abnormal  ST  segment  depression 
on  symptom-limited  stress  testing  have  no  coronary  artery  disease  found  on 
cardiac  catheterization  with  left  ventriculography  and  selective  coronary 
angiography.  We  have  also  found  advanced  coronary  artery  disease  in  three 
dozen  asymptomatic  flyers  with  normal  or  only  borderline  abnormal  ST 
depression  on  their  treadmill  test. 

In  a  series  of  four  articles  published  in  Circulation  in  1973  and  1974, 
Fred  Kornreich,  M.D.,  of  the  Free  University  of  Brussels,  pointed  out  that 
neither  the  Frank  lead  vectorcardiogram  nor  the  standard  12-lead  electrocar¬ 
diogram  contain  all  the  information  projected  by  the  heart  to  the  body  sur¬ 
face.  Dr.  Kornreich  found  that  with  24  leads  and  a  single  transformation 
formula  for  all  body  types,  one  could  closely  reconstruct  a  complete  body 
surface  map.  Subsequently,  Robert  Lux,  Ph.D.,  at  the  University  of  Utah, 
Roger  Barr,  Ph.D.,  at  Duke  University,  and  others  have  shown  that  at  least  24 


electrodes  are  needed  to  capture  all  the  information  projected  by  the  heart 
to  the  body  surface.  What  remains  to  be  determined  is  whether  24  electrodes 
are  necessary  to  detect  all  disqualifying  asymptomatic  coronary  artery 
disease  in  the  USAF  rated  force.  In  1980,  Dr.  Kornreich  reported  at  the 
Engineering  Foundation  Conference  on  Computerized  Interpretations  of  the  ECG 
that  adding  just  four  electrodes  to  the  standard  12-lead  ECG  for  a  total  of 
14  electrodes  significantly  improved  the  sensitivity  and  specificity  of 
bonafide  myocardial  infarctions  over  just  the  standard  12-lead  ECG  or  just 
the  Frank  lead  VCG.  At  the  same  meeting,  Frank  Yanowitz,  M.D.,  of  the 
University  of  Utah,  reported  that  using  32  leads  doubled  the  sensitivity  of 
stress  tests  for  cardiac  catheterization-proven  significant  coronary  artery 
disease  compared  to  using  the  standard  12-lead  system  for  stress  testing.  He 
also  suggested  that  not  all  32  leads  were  necessary,  but  needed  more  data  to 
prove  this.  In  1980,  at  the  Computers  in  Cardiology  Conference,  R.  von  Essen 
of  the  Helmholtz  Institute  for  Biomedical  Engineering  in  West  Germany, 
reported  that  a  48-precordial-lead  stress  test  had  a  96%  sensitivity  and  95% 
specificity  compared  with  56%  sensitivity  and  95%  specificity  of  the  12-lead 
system  for  stress  testing.  Other  investigators  have  reported  similar 
improvements  in  sensitivity  while  maintaining  95%  specificity  by  doing 
surface  potential  mapping  with  stress  testing. 


8 


NOMENCLATURE  OF  REGIONAL  LEFT  VENTRICULAR  SEGMENTS 


Considerable  ambiguity  exists  in  the  current  literature  about  the  nomen¬ 
clature  and  definition  of  the  regional  or  segmental  anatomy  of  the  heart. 
Throughout  this  document  we  will  be  using  the  left  ventricular  (LV)  segmental 
subdivision  nomenclature  shown  in  Figure  1.  This  subdivision  is  a  modifica¬ 
tion  of  the  Ideker  subdivision  reported  in  a  number  of  recent  papers  (18-23), 
in  which  the  octant  circumferential  subdivisions  of  the  left  ventricle  along 
its  long  axis  are  combined  by  twos  into  quadrants.  The  LV  is  further  sub¬ 
divided  from  apex  to  base  by  passing  two  planes  through  it  at  right  angles  to 
the  long  axis  of  the  LV  in  such  a  way  as  to  divide  the  internal  long  axis 
into  equal  thirds.  This  system  produces  a  12-segment  subdivision  of  4  walls 
(anterior,  superior,  posterior,  and  inferior)  with  3  subdivisions  each 
(apical,  mesial,  and  basal)  as  shown  in  Figure  1. 


Figure  1,  Twelve-segment  LV  subdivision  and  nomenclature  used  in  this 

report.  The  advantage  of  using  this  specific  subdivision  is  that 
it  divides  the  left  ventricle  into  segments  that  are  in  close 
approximation  to  the  usual  blood  supply  to  each  of  the  regions, 
i.e.,  segments  1-7,  and  10  by  the  left  anterior  descending 
coronary,  8  and  9  by  the  right,  and  11  and  12  by  the  circumflex. 


SOME  WORKING  DEFINITIONS 


Activation,  excitation,  and  depolarization  of  the  ventricle  (which 
produces  the  QRS  of  the  ECG)  will  be  used,  throughout  this  report, 
interchangeably.  The  same  applies  to  recovery  and  repolarization,  the  ST-T 
portion  of  the  ECG.  End  diastole  will  be  used  for  that  time  in  the  cardiac 
cycle  when  maximal  diastolic  relaxation  has  occurred  just  prior  to  the  onset 
of  ventricular  contraction.  The  QRS  of  the  ECG  is  recorded  at  this  time. 
End  systole  will  be  taken  to  mean  that  time  when  the  ventricular  contraction 
has  reached  its  maximum,  the  slope  of  the  ventricular  pressure  curve  is  0, 
and  ventricular  relaxation  will  soon  commence,  followed  by  aortic  and  pul¬ 
monic  valve  closure,  and  the  return  of  the  ventricular  pressure  curve  to 
diastolic  levels.  Mid  systole  is  that  time  midway  between  end  diastole  and 
end  systole. 
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BACKGROUND  AND  PREVIOUS  WORK 


Ventricular  Depolarization 


Extensive  work  has  been  done  in  observing  and  measuring  the  complex  bio¬ 
electric  phenomenon  of  depolarization  of  myocardial  cells.  Over  a  century 
ago,  Burdon-Sanderson  and  Page  used  mercury  capillary  electrometers  to  mea¬ 
sure  in  some  detail  the  electric  current  generated  by  both  auricles  and 
ventricles  (54,55).  Giving  credit  to  previous  work  in  the  mid  1800‘s  by 
Kolliker  and  Muller  (56),  Cyon  (57),  Donders  (58),  Engellman  (59),  and 
Marchand  (60),  they  observed  that  the  positivity  of  the  electrical  deflection 
that  preceded  contraction  of  the  frog  (and  turtle)  heart  was  related  to  the 
location  of  the  applied  electrodes  as  related  to  the  long  axis  of  the  ventri¬ 
cles.  They  recorded  a  terminal  wave,  later  labeled  the  T  wave  by  Einthoven, 
and  noted  that  it  coincided  with  ventricular  relaxation.  They  found  that 
warming  or  cooling  changed  the  timing  and  the  configuration  of  this  wave  and 
that  an  electrical  stimulus  delivered  before  the  terminal  wave  did  not  pro¬ 
duce  a  second  primary  electrical  wave  and  contraction,  whereas  a  stimulus 
delivered  after  the  terminal  wave,  but  before  the  next  primary  wave,  did. 

The  invention  of  the  string  galvanometer  at  the  turn  of  the  century  and 
its  development  by  Einthoven  into  a  routine  experimental  tool  (1,61,62)  led 
to  a  great  deal  of  experimental  work  in  animal  models.  Notable  was  the 
careful  work  in  the  dog  by  Lewis  (2-4)  who  described  normal  and  abnormal 
activation  sequences  based  on  epicardial  measurements  with  the  string  gal¬ 
vanometer.  Lewis1  version  of  the  sequence  of  normal  ventricular  depolariza¬ 
tion  and  its  extrapolation  to  man  was  the  basic  model  used  to  interpret 
clinical  ECGs  for  the  next  30  years.  A  number  of  workers  explored  the  role 
of  the  Purkinje  conduction  system  in  controlling  the  activation  sequence 
(63-69).  Conflicting  data  was  generated  which  led  to  a  misnaming  of  right 
and  left  bundle  branch  block  until  a  review  paper  by  Wilson  (70)  that 
included  convincing  new  work  clarified  the  confusion.  The  development  of  the 
precordial  leads,  first  the  CF  leads,  and  later  the  V  leads  using  the  central 
terminal  of  Wilson  (71-73)  added  a  much  needed  horizontal  plane  dimension  to 
the  ECG.  This  development  helped  clarify  the  localization  of  bundle  branch 
blocks  and  local  infarction  as  well.  The  development  of  various  vectorcar- 
diographic  lead  systems  soon  followed,  which  further  refined  the  early 
Einthoven  hypothesis  that  the  heart  could  be  well  represented  by  a  single 
fixed-locus  equivalent  dipole  (74-78). 

Papers  by  Prinzmetal  and  coworkers  (79,80)  and  Sodi-Pal lares  and  cowor¬ 
kers  (81-84)  recording  with  intracavitary  and  plunge  electrodes  in  the  dog 
heart  added  information  about  the  electrical  behavior  of  the  heart  with 
ischemia  and  infarction.  However,  the  sequence  of  intramural  excitation  in 
the  presence  of  infarction,  or  even  in  normal  activation,  remained  sketchy 
due  to  the  limited  number  of  electrode  sites  on  these  electrodes.  The 
detailed  intramural  mapping  of  the  wave  of  excitation  of  the  mammalian  heart 
using  small  diameter,  multi-point  needle  electrodes  was  reported  in  a  number 
of  papers  by  Durrer  and  coworkers  (9-13)  and  in  the  classic  papers  by  Scher 
and  his  coworkers  (5-8).  Soon  thereafter,  Durrer,  Van  Dam,  and  associates 
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presented  the  details  of  human  activation  sequence  using  revitalized  human 
hearts  of  accident  victims  (85).  These  latter  papers  have  become  the  basis 
of  a  number  of  simulations  of  ventricular  depolarization  in  man,  including 
those  of  Okajima  (86,87),  Horachek  et  al.  (88),  and  Arntzenius  (89),  as  well 
as  our  own. 

In  the  early  1960‘s,  we  began  to  explore  the  possible  role  of  computer 
simulations  of  the  electric  field  of  the  heart  in  man.  Our  first  effort  was 
to  simulate  the  human  vectorcardiogram  (90).  We  assumed  the  heart  to  be  an 
equivalent  dipole  current  source  and  the  human  activation  sequence  to  be 
similar  to  that  demonstrated  by  Scher  in  dogs.  We  also  assumed  a  uniform 
current  generator  surface  in  the  local  region  of  a  multiple  segment  heart. 
We  gave  each  segment  a  dipole  moment  time  history  based  on  these  assumptions 
and  the  Scher  activation  data,  and  a  direction  normal  to  the  average  wave 
fronts  in  the  local  region.  When  normal  activation  was  simulated,  normal 
vectorcardiograms  resulted.  When  dipoles  representing  the  left  or  right 
ventricle  were  increased  in  dipole  moment,  as  would  be  expected  with  hyper¬ 
trophy  of  either  ventricle,  VCGs  typical  of  these  abnormalities  were  pro¬ 
duced,  and  when  one  dipole  or  more  were  removed  to  simulate  the  loss  of  8 %  or 
more  of  the  left  ventricle,  VCG  loops  (and  simulated  12-lead  ECGs)  were 
produced  that  were  quite  typical  of  the  focal  infarct  simulated.  Moreover, 
the  larger  the  infarct  simulated,  the  greater  the  VCG/ECG  change.  When 
Durrer's  human  activation  data  became  available  at  about  this  time,  we  found 
it  very  similar  to  our  simulation.  We  upgraded  our  model  with  this  measured 
data,  and  with  distance  and  boundary  effects  (15,16).  A  year  later,  we 
added  the  Gelernter-Swihart  simulation  of  the  inhomogeneous  human  male  torso 
including  intracavi tary  blood  mass  (91).  Finally,  we  added  a  1  mm3  grid 


SAGTTTAL-  <V«rt) 


Isochronous  time  lines  on  a  sagittal  cross-section 
corresponding  to  the  same  propagation  of  fig.  3. 


LEVEL  -30  mm  (Hors.) 

Isochronous  time  lines  every  10  milliseconds  represent 
the  position  of  the  wavefront  on  a  horizontal  cross-section  of 
the  cardiac  geometry.  Purkmie  velocity  is  four  limes  faster 
than  cell-io-cell  conduction  velocity  of  40  cm/sec. 


Figure  2.  A  transverse  view  along  the  long  axis  of  the  left  ventricle  is 

shown  to  the  left,  and  a  short  axis  view,  corresponding  to  the  30- 
mm  cross  section  identified  on  the  long  axis  view, is  shown  to  the 
right.  Isochronous  time  lines  every  10  milliseconds  are  shown  on 
both  views. 
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propagation  model  of  excitation  (Figure  2)  which  used  the  actual  digitized 
geometry  of  a  normal  male  human  heart  and  included  a  simulation  of  Purkinje 
conduction  (92,93).  We  now  had  in  the  simulation  all  the  first-order  vari¬ 
ables  known  to  influence  the  human  ECG,  and  for  all  intents  and  purposes,  the 
computer  simulation  behaved  like  the  electric  field  of  a  normal  human  male 
subject.  The  parameters  that  influence  wave  shape,  such  as  torso  geometry, 
resistivities,  heart  location,  and  activation  sequence  of  normal  (Figure  3) 
or  abnormal  (Figure  4)  conduction  with  or  without  hypertrophy  or  infarction, 
could  be  varied  at  will  within  physiologically  and  anatomically  reasonable 
limits.  The  complex  interactions  of  these  variables  could  be  examined  in  a 
situation  where  all  the  parameters  used  and  the  assumptions  made  are 
precisely  known. 

In  the  early  stages  of  the  experimentation  with  the  simulations  that 
could  mimic  all  the  classic  infarcts,  it  became  clear  why  conventional  ECG 
criteria  are  insensitive  to  many  infarcts.  In  a  number  of  autopsy  series 
(94-96),  the  ECG  was  found  predictive  of  the  infarcts  found  in  only  about  55% 
to  60%  of  those  confluent  scars  that  were  2  cm  or  greater  in  size  at  post 
mortem  examination.  Conventional  ECG  criteria  for  infarction  involve  the 
finding  of  abnormal  Q  waves  or  abnormal  loss  of  initial  R  waves  especially  in 
precordial  leads.  In  the  simulation,  small  infarcts  in  the  typical  locations 
produced  definite  changes  when  compared  to  pre-infarct  ECG/VCGs,  but  often 
failed  to  produce  the  classical  changes.  Larger  infarcts  in  these  same 
locations  generally  did  produce  typical  changes.  On  the  other  hand,  about 
25%  of  infarcts  found  at  autopsy  involve  the  base  of  the  left  ventricle  which 
is  depolarized  after  40  msec.  Simulated  infarction  of  these  regions  produced 
mid  to  late  QRS  changes  and  no  initial  changes  at  all  that  would  be  the 
typical  or  classic  ECG  changes  of  infarct.  We  designed  one  of  our  early 
experiments  (97)  to  explore  the  contention  by  Sodi-Pal lares  (84,98)  and 
others  (79)  that  large  portions  of  the  basal  septum  and  free  wall,  as  well  as 
most  of  the  subendocardium,  is  electrically  silent.  We  found  that  neither 
extensive  experimental  work  with  the  dog  (99)  and  the  baboon,  nor  extensive 
work  with  the  model  that  contained  the  details  of  electrophysiology  and 
anatomy  (14-16)  supported  this  notion.  In  fact,  while  some  regions  deep  in 
the  chest  showed  an  expected  lower  amplitude  of  change  in  surface  ECG/VCGs 
and  surface  ECG  maps,  no  region  was  electrically  silent.  To  the  contrary, 
simulated  infarcts  in  all  regions  of  the  heart  produced  clear-cut  changes 
when  compared  to  prior  ECG  data.  Moreover,  the  degree  of  change  was 
proportional  to  the  extent  of  infarct  put  into  the  model. 

We  next  explored  the  question  of  what  size  of  infarct  in  any  of  the 
various  regions  of  the  heart  would  produce  detectable  changes  in  optimally 
recorded  12-lead  ECGs  and  VCGs  by  both  a  corrected  [McFee-Parungao  (100)]  and 
an  uncorrected  [modified  Grishman  Cube  (77)]  lead  system.  When  incorporated 
into  the  propagation  simulation  and  the  full  torso  model  (101),  lesions  7  mm 
to  10  mm  in  diameter  and  3  mm  thick  produced  serial  notches  and  deformities 
in  the  QRS  well  above  the  noise  in  well-recorded,  signal-processed  ECG  data 
(15).  Lesions  15-20  mm  in  diameter  and  4-5  mm  thick  produced  quite  unequivo¬ 
cal  changes  (16).  Examples  of  these  changes  are  shown  in  Figures  5  and  6  and 
are  shown  in  more  detail  in  refs.  15,  16,  22,  and  101. 
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VENTRICULAR  PROPAGATION  MODEL  WITH  NORMAL  CONDUCTION 


Figure  3 


A  moderate  to  large  anterioapical  infarct  that  involves  24%  of  the 
left  ventricle  is  shown  in  the  diagonal  lined  shaded  areas  in  the 
apical  and  mid-ventricular  cross  sections.  It  is  incorporated 
into  the  model  along  with  the  normal  activation  as  shown  in  Fig. 
2.  The  early  activation  on  the  mid  left  septal  surface  and  up 
onto  the  endocardial  surface  of  the  superior  and  inferior  walls  is 
shown  in  solid  black  in  this  diagram.  Isochronous  time  lines  at 
10-msec  intervals  are  shown  as  the  excitation  wave  propagates 
outward  from  the  normal  early  start  points  on  the  endocardium.  The 
40-msec  isochrone  is  shown  in  darker  shaded  line  labeled  40  and 
best  seen  in  the  uninfarcted  basal  section.  As  shown  here  it  is 
common  for  infarcts  in  the  human  heart  to  involve  the  endocardium 
more  extensively  than  the  epicardium  and  to  leave  a  thin  "peel"  of 
relatively  normal  myocardium  in  the  epicardium  over  the  infarct. 
When  this  happens,  the  excitation  reaches  the  epicardium,  as  in 
this  case,  from  the  endocardium  around  the  infarct.  Except  for 
this  small  area  of  "peri-infarction"  delay,  the  latest  area  acti¬ 
vated  was  the  posterior  wall  up  near  the  AV  groove,  and  the  out¬ 
flow  tract  of  the  right  ventricle,  not  shown  in  any  of  these  three 
cross  sections.  For  more  specific  details  about  how  various  sized 
infarcts  that  were  simulated  in  each  coronary  artery  distribution 
affected  the  12-lead  ECG  and  the  McFee  VCGs  on  an  inhomogeneous 
male  torso,  see  the  text,  and  particularly  ref.  22. 
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SUPERIOR  AND  INFERIOR  FASCICULAR  BLOCK  IN  THE  PROPAGATION  MODEL 


Figure  4.  Superior  fascicular  block  was  simulated  in  this  experiment  by 
omitting  the  start  points  in  the  superior  wall  of  this  heart. 
Activation  reached  this  region  as  shown  here  after  about  a  20  msec 
delay,  arriving  there  via  an  intact  peripheral  endocardial 
Purkinje  network.  Such  an  excitation  produces  20-30  msec  increase 
in  total  QRS  duration,  a  shift  in  early  forces  inferior  (producing 
significant  Q's  in  I  and  AVL),  and  late  forces  superior.  Because 
of  major  unopposed  late  activation  fronts  superior  (and  poster¬ 
ior),  there  is  a  shift  of  the  frontal  mean  axis  superior  to  above 
-40  degrees.  Inferior  fascicular  block  in  this  experiment  pro¬ 
duced  similar  degree  of  total  QRS  prolongation,  initial  vectors 
superior  for  30  msec  and  a  prominent  unopposed  late  front  inferior 
and  posterior.  The  net  effect  of  these  changes  was  to  shift  mean 
frontal  axis  inferior  to  +70  degrees.  It  took  an  unusual  variant 
of  normal  excitation  to  produce  "right  axis"  (i.e.,  beyond  90  to 
100  degrees  in  the  frontal  plane  with  a  block  in  the  inferior 
fascicle).  We  believe  that  this  accounts  for  the  relative  rarity 
of  this  diagnosis  using  conventional  criteria  for  its  diagnosis  by 
ECG/VCG.  Specifically,  the  criteria  for  the  diagnosis  of  inferior 
fascicular  block  appear  to  be  too  stringent.  In  a  critical  review 
of  our  data  base,  we  found  that  this  conduction  abnormality  is 
about  as  common  as  superior  fascicular  block  (22)  and  is  seen  most 
commonly  in  association  with  inferior  infarction  from  right 
coronary  occlusion. 


Figure  5.  A  small  area  of  infarction  in  the  outer  half  of  segment  11  is 
simulated  here.  This  simulation  was  done  by  placing  a  dipole  with 
a  negative  directed  function,  represented  by  the  shaded  area 
above,  at  the  same  coordinate  location  and  given  the  direction 
cosines  as  segment  11.  The  computation  then  sees  both  dipoles  as 
one  and  the  resultant  will  be  the  double  peaked  curve  shown.  The 
infarcted  area  shaded  above  is  about  10%  of  the  area  under  the 
curve  representing  the  dipole  moment  history  of  segment  11.  Thus, 
this  vectorcardiographic  change  was  produced  by  a  simulated 
infarct  representing  about  2  g  of  myocardium. 


The  validity  of  these  findings  was  tested  in  a  series  of  patients  with 
known  coronary  artery  disease  in  whom  serial  coronary  and  biplane  ventricular 
angiograms  were  done  as  a  part  of  a  long-term  follow-up  study.  All  had 
serial  12-lead  ECGs,  and  Cube  and  McFee  VCGs  digitized  at  a  1  kHz  sampling 
rate,  cycle  selected  in  expiration  time  aligned,  and  averaged  to  reduce  noise 
to  the  3  microvolt  peak-to-peak  range.  The  presence  or  absence  of  signifi¬ 
cant  serial  change  as  read  by  two  blinded  readers  was  predictive  with  a  98% 
accuracy  of  the  presence  or  absence  of  significant  serial  ventricul ographic 
change  as  read  by  2  different  blinded  readers  (17).  These  findings  along 
with  the  simulations  led  to  the  development  of  a  54  criteria-32  point  ECG 
infarct  scoring  system  where  each  point  represents  3%  infarction  of  the  left 
ventricle.  The  point  scoring  system,  shown  in  detail  in  Table  1,  has  been 
extensively  validated  in  a  series  of  autopsy  studies  (18-20,22),  and  by 
clinical  nuclear  and  contrast  angiography  studies  by  Wagner,  Ideker,  and 
associates  (21,102).  While  the  usual  12-lead  ECGs  used  for  these  studies 
were  not  optimal  by  the  criteria  of  high-fidelity,  low-noise  records  as 
defined  above,  the  good  correlations  achieved  (see  Figure  7)  indicate  that 
even  in  these  poor  quality  records,  clinically  relevant  quantitative 
information  about  infarct  size  and  residual  ventricular  function  is  present. 
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Figure  6. 


A  tiny  infarct  representing  less  than  a  6-mm  lesion  in  segment  8 
at  the  apex  of  the  left  ventricle  is  shown  here.  When  a  control 
record  was  available  for  comparison,  a  simulated  small  infarct  of 
this  size  produced  recognizable  changes  in  any  of  the  20  dipole 
locations.  These  data  cast  doubt  on  the  concept  of  electrically 
"silent  area"  of  myocardium,  and  suggest  instead  that  what  is 
needed  is  to  refine  our  diagnostic  criteria  and  our  measuring 
tools  so  that  they  are  sensitive  to  this  degree  of  change  in  the 
human  cardiac  generator. 


Clinical ly  "Silent"  Infarcts  in  the  Air  Force  Pilot 


Regarding  the  case  finding  of  clinically  "silent"  infarcts  in  the  asymp¬ 
tomatic  Air  Force  pilot  population,  these  data  indicate  that  in  this  popula¬ 
tion  where  serial  ECGs  are  routine  and  high-fidelity  VCGs  are  available  on  a 
large  number,  significant  infarction  (i.e.,  infarcts  greater  than  1.5  cm  in 
diameter)  should  be  routinely  detectable.  Since  some  20-30%  of  infarcts  in  a 
civilian  population  are  clinically  silent,  and  since  this  percentage  in  the 
pilot  population  can  be  expected  to  be  about  twice  as  high,  it  follows  that 
approximately  50  infarcts  occur  each  year  in  rated  pilots  and  navigators  that 
are  undetected  clinically.  Most  of  these  infarcts  are  probably  small  to  mode¬ 
rate  in  size  since  most  infarcts  are  in  this  size  range.  Thus,  well  over 
half  of  these  would  not  have  classic  ECG  changes  on  a  routine  serial  ECG 
screening.  With  careful  attention  to  recording  serial  high-fidelity  ECG/VCG 
data  and  to  criteria  involving  the  entire  QRS,  including  significant  deform¬ 
ities  not  associated  with  classic  Q  wave  changes,  it  is  to  be  expected  that 
fully  90%  or  more  of  these  infarcts  can  be  found.  Since  a  small  number  of 
false  positives  are  inevitable  from  any  such  review  of  data,  we  propose  that 
the  serial  tracings  be  screened  for  major  significant  non-Q  wave  changes, 
indicative  of  at  least  a  moderate-si  zed  new  infarct  by  the  ECG  Point  score 
(see  Table  1)  or  by  serial  VCG  change.  Moderate  infarcts  would  be  5-7  ECG 
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Figure 


PANEL  A 


PANEL  B 


PATHOLOGIC  INFARCT  SIZE  AND  ECG  CHANGES 
(ALL  LOCATIONS) 


7.  Panel  A  plots  new  infarct  size  on  serial  comparisons  of  VCG 
change  (or  lack  of  it)  and  contrast  Vgrams  in  77  consecutive 
serial  coronary  angiographic  studies  (17).  The  upper  line  is  the 
regression  line,  and  the  lower  one  is  line  of  unity.  Panel  B 
shows  the  correlation  between  our  ECG  point  score  for  infarct  size 
and  ejection  fraction  (E.F.)  in  patients  with  chronic  angina 
and/or  infarction  (22).  The  2  extra  lines  represent  1SD  on  either 
side  of  the  regression  line.  Panel  C  relates  ECG  points  and 
nuclear  E.F.  in  43  fatients  3  weeks  after  an  acute  infarction,  the 
dotted  lines  =  2SD;  Palmeri  et  al.  (102).  Panel  D  plots  ECG  point 
score  and  infarct  size  at  quantitative  planometric  pathology  in 
the  1st  53  patients  reported  by  Ideker  et  al.  (19-23). 


19 


points,  or  15-20%  of  the  left  ventricle.  This  would  include  a  similar  degree 
of  VCG  change  by  published  criteria.  The  intention  is  to  approach  a  100% 
likelihood  that  pilots  with  these  changes  on  their  ECGs  would  have  signifi¬ 
cant  coronary  artery  disease  on  coronary  angiograms. 


Cardiac  Intracellular  Action  Potentials 

The  development  of  the  fine  glass  capillary  microelectrode  by  Ling  and 
Gerard  in  1949  (103)  led  to  a  new  era  in  the  exploration  of  intracellular 
action  potentials  from  electrically  active  cells.  The  small  (1  ym  or  less) 
tip  size  allowed  for  micropuncture  of  a  single  cell  without  significantly 
changing  its  membrane  electrophysiology.  In  1951,  Woodbury,  Hecht  and 
Christopherson  recorded  single  ventricular  fiber  action  potentials  from  the 
frog  using  this  electrode  (104);  with  Brady,  this  group  examined  the  effects 
of  cardiac  glycosides  and  other  interventions  on  this  frog  ventricular  fiber 
preparation,  and  later  measured  transmembrane  action  potentials  from  the 
human  heart  (105-107).  Seminal  work  relating  the  ventricular  transmembrane 
action  potential  duration  and  wave  shapes  to  refractory  period  and  the  T  wave 
of  the  ECG  was  done  in  the  1950 * s  by  Hoffman,  Kao,  and  Suckling  (108,109). 
This  work  culminated  in  the  publication  in  1960  of  the  classic  monograph  by 
Hoffman  and  Cranefield,  "Electrophysiology  of  the  Heart"  (110).  This  publi¬ 
cation  extensively  reviewed  transmembrane  action  potentials  from  various 
cardiac  tissues  and  the  effect  of  various  interventions  using  cardioactive 
drugs,  electrolytes,  hypoxia,  and  injury  on  these  Cardiac  ce^s. 

The  shape  of  the  action  potiental,  however,  can  be  recorded  relatively 
free  of  the  problems  of  muscle  movement  in  the  intact  heart  if  external 
suction  electrodes  are  used  (111).  Six  years  before  the  work  of  Ling  and 
Gerard  with  the  cell  puncture  with  microelectrodes,  Schafer  et  al.  had  used 
the  suction  electrode  to  study  the  T  wave  (112).  In  1959,  Hoffman  and 
Cranefield,  working  with  Lepeschkin,  Surawicz,  and  Herrlich,  demonstrated 
that  the  waveform  of  the  action  potential  by  the  two  methods  corresponded 
closely  although,  as  might  be  expected,  the  magnitude  of  the  potential 
recorded  with  the  suction  electrode  was  smaller  (113).  The  main  disadvantage 
of  the  suction  electrode  is  that  it  must  be  applied  to  the  surface  of  the 
intact  heart  and  therefore,  cannot  be  used  to  gather  information  about  gra¬ 
dients  of  action  potential  durations  and  sequence  of  recovery  in  the  intra¬ 
mural  layers  of  the  heart.  Another  limitation  of  the  method  is  that  it  must 
also  be  used  with  the  chest  open,  and  repolarization  is  known  to  be  quite 
sensitive  to  temperature  changes.  Using  this  method  and  attempting  to  mini¬ 
mize  the  temperature  effects,  Autenrieth,  Surawicz,  and  Kuo  were  able  to 
demonstrate  systematic  differences  in  the  duration  of  action  potentials  on 
the  surface  of  the  left  ventricle  of  the  dog  (114).  The  longest  duration  was 
at  the  apex,  and  the  shortest  was  at  the  posterior  base.  They  also  noted 
that  in  most  cases,  the  surface  action  potentials  ended  well  before  the  end 
of  the  T  wave,  implying  that  some  tissue  deeper  in  the  wall  of  the  heart,  not 
accessible  to  the  electrode,  repolarized  later. 

Based  on  the  findings  of  Hoffman  and  Cranefield  that  refractory  periods 
measured  by  classic  methods  corresponded  closely  to  the  end  of  the  repolari¬ 
zation  phase  of  the  intracellular  action  potentials,  van  Dam  and  Durrer 
mapped  functional  refractory  periods  with  intramural  electrodes  in  the  open¬ 
chested  dog  (115,116).  They  found  that  the  endocardium,  while  depolarizing 


first,  recovered  refractoriness  latest;  mid  wall  intramural  areas  recovered 
earliest;  and  the  epicardium  is  intermediate.  Burgess,  Abildskov,  Green, 
Lux,  Millar,  Vincent,  Wyatt,  Yanowitz  and  others  of  the  Utah  group  published 
an  important  series  of  papers  in  the  last  15  years  repeating  this  work  in  the 
closed-chested  dog  and  concluded  that  there  is  likely  a  linear  gradient  of 
ventricular  recovery  times  from  epicardium  to  endocardium  with  the  epicardium 
repolarizing  earliest  (117-132).  They  found  that  the  interventricular  septum 
behaves  like  the  free  wall  of  the  left  ventricle.  They  also  established  that 
there  is  a  gradient  in  functional  refractory  periods  from  apex  to  base  with 
the  apical  regions  having  the  longer  refractory  periods. 

Finally,  extracellular  electrode  recordings,  both  bipolar  and  unipolar, 
have  been  used  to  map  excitation  in  detail.  Their  advantage  over  even  multi¬ 
ple  intracellular  recordings  is  that  they  map  the  interactions  of  the  whole 
system  of  depolarizing  and  repolarizing  cells.  Spach  and  his  coworkers  made 
unipolar  potential  maps,  not  only  with  electrodes  on  the  torso  surface,  but 
with  epicardial  and  intramural  electrodes  in  animals  allowed  to  recover  for 
several  days  before  the  measurements  were  made  (133-135).  From  some  300 
electrodes  implanted  in  each  of  several  dogs,  they  were  able  to  reconstruct 
detailed  maps  of  depolarization  that  agreed  with  earlier  work  by  them  and  a 
number  of  others  mentioned  in  this  report.  In  addition,  they  were  able  to 
map  the  recovery  process  in  considerable  detail.  This  preparation  has  the 
advantage  that  recovery  is  mapped  in  the  same  animal  with  typically  normal 
epicardial  and  body  surface  ECG  potential  distributions  for  both  depolariza¬ 
tion  and  repolarization.  Careful  review  of  these  papers  confirms  in  detail 
the  refractory  period  data  of  Burgess  et  al.,  in  that  the  earliest  recovery 
potentials  were  recorded  in  the  right  ventricular  side  of  the  interventricu¬ 
lar  septum,  with  epicardial  basal  regions  beginning  to  show  recovery  poten¬ 
tials  soon  thereafter.  They  found  that  the  endocardium  was  positive  with 
reference  to  the  epicardium  throughout  recovery  with  the  endocardial  regions 
at  the  apex  being  the  last  area  to  show  significant  recovery  potentials. 
Spach  comments  that  this  data  is  consistent  with  the  linear  gradient  in  epi- 
to  endocardial  recovery  times  postulated  by  Burgess,  Abildskov  and  coworkers, 
and  is  inconsistent  with  the  proposal  by  van  Dam  and  Durrer  that  the  mid  wall 
of  the  left  ventricle  recovered  before  either  the  epicardium  or  the  endocar¬ 
dium.  These  data  for  depolarization  and  repolarization  are  shown  for  our  12 
segment  left  ventricular  subdivision  in  Figures  8  and  9. 


Modelling  ot  Ventricular  Kepolarization 

A  comprehensive  critical  review  entitled  "The  Origin  of  the  T-Wave"  was 
published  in  CRC  Critical  Reviews  in  Bioengineering  in  November  1980  by 
Kootsey  and  Johnson  with  M.  Spach  as  referee  (136).  This  monograph  will  be 
excerpted  where  directly  applicable,  but  no  attempt  will  be  made  to  replicate 
the  entire  contents.  Since  this  monograph  was  published,  little  new  has  been 
added  to  the  literature  that  relates  to  the  topic  of  simulation  of  recovery 
and  the  reader  is  urged  to  consider  the  Kootsey  and  Johnson  paper  as  an 
integral  part  of  this  technical  report. 

While  all  models  of  repolarization  are  indebted  to  some  extent  to  the 
notion  of  Ventricular  Gradient  summarized  in  the  introduction  and  to  the 
relating  of  action  potential  durations  and  functional  refractory  periods 
described  in  the  previous  section,  the  first  attempt  to  describe  a 


figure  8.  Activation  as  shown  in  the  revitalized  human  hearts  by  Durrer  and 
van  Dam  and  by  us  in  the  chimpanzee  and  baboon  with  detailed 
intramural  multipoint  mapping  is  depicted  here.  Colored  iso¬ 
chrones  for  each  10  msec  are  those  defined  by  the  expanded  color 
code  at  the  lower  right.  They  depict  the  normal  activation  incor¬ 
porated  in  our  fine  grid  propagation  model  and  the  12-segment 
subdivision  used  in  the  development  of  54  criteria/32  point  ECG 
infarct  size  scoring  system  shown  in  Table  1. 


n.jur-p  J.  The  sequence  of  ventricular  recovery  is  depicted  here  adapted  from 
the  functional  refractory  period  data  of  Burgess,  Abildskov,  and 
coworkers,  and  from  the  detailed  intramural  mapping  of  this  phase 
of  the  cardiac  cycle  by  Spach  and  coworkers  in  the  chronic  closed¬ 
chested  dog  described.  Recovery  occurs  much  more  slowly  than 
excitation  and  begins  soon  after  depolarization,  as  indicated  by 
the  action  potentials  shown  in  the  lower  right.  The  color  coding 
indicates  the  timing  of  the  peak  of  the  slope  of  the  rapid  recov¬ 
ery  phase  in  each  region.  According  to  Spach,  as  shown  here,  the 
earliest  recovery  potentials  appear  in  the  right  side  of  the 
septum,  recovery  proceeds  more  rapidly  at  the  base  and  from  epi- 
cardium  to  endocardium  throughout.  The  latest  recovery  potentials 
are  seen  in  the  endocardium  near  the  apex. 


quantitative  model  of  the  T  wave  was  in  1966  by  Harumi,  Burgess,  and 
Abildskov  (50).  This  early  version  was  a  2-dimensional  one  in  which  the 
cardiac  sources  moved  through  a  simple  ventricular  cross  section.  The  local 
repolarization  was  assumed  to  be  driven  by  the  activation  sequence  in  that 
cross  section.  A  transmural  gradient  of  action  potential  duration  was 
assumed  (as  shown  in  Figure  8)  with  the  epicardial  durations  being  shorter 
than  the  endocardial.  The  magnitude  of  this  duration  gradient  was  adjusted 
so  that  the  T  vector  loops  computed  by  the  solid  angle  theorem  conformed  to 
the  experimentally  recorded  ones.  The  group  at  the  University  of  Utah, 
including  Abildskov,  Burgess,  Green,  Lux,  and  others  (117-132),  has  been 
active  for  the  last  several  years  in  the  investigation  of  the  phenomenon  of 
ventricular  recovery.  They  have  established  themselves  at  the  forefront  of 
this  activity  worldwide.  While  they  have  been  very  active  investigating 
regional  recovery  with  cooling,  coronary  ligation,  focal  ischemia,  and  ven¬ 
tricular  pacing,  in  general,  they  have  limited  their  modelling  efforts  to  the 
simple  2D  model . 

Based  on  the  data  from  the  University  of  Utah  Group,  Toyama,  Okajima,  and 
associates  (137),  and  Hori,  Inoue,  and  their  coworkers  (138,139,140),  have 
reported  various  experimental  simulations  using  strand,  concave  rectangular 
solids,  and  spherical  models  to  explain  repolarization.  Little  attempt  to 
model  the  volume  conductor  was  made.  The  numerical  experiments  concerned 
solid  angle  views  of  the  surface  in  an  infinite  homogeneous  medium  and  exam¬ 
ined  the  degree  of  gradient  in  action  potential  durations  needed  to  simulate 
a  normal  T  vector.  They  established  in  these  simulations  that  the  gradients 
from  epicardium  to  endocardium  that  produced  upright  T  waves  conformed  to 
those  observed  in  the  animal  laboratory  by  Abildskov,  Burgess,  and  coworkers. 

Thiry  and  Rosenberg  (51,141,142)  were  the  first  to  attempt  a  3- 
dimensional  model  in  the  form  of  the  cardiac  geometry  to  model  both  excita¬ 
tion  and  the  T  wave.  Eleven  current  dipoles  were  fixed  in  location  and 
direction  to  represent  each  of  11  regions  of  the  heart.  A  model  of  individ¬ 
ual  action  potential  wave  shapes  was  developed  where  the  duration,  amplitude, 
and  downslope  could  be  individually  varied  to  simulate  the  Burgess-Abildskov 
T  model  and  their  refractory  period  duration  distributions  in  the  wall  of  the 
left  ventricle.  As  a  matter  of  mathematical  convenience,  each  local  myocar¬ 
dial  region  was  represented  by  a  sphere.  With  this  formulation,  equations 
could  be  written  such  that  different  action  potential  waveshapes  would  repre¬ 
sent  endocardial,  mural  myocardial,  and  epicardial  regions  within  this  spher¬ 
ical  "cell".  In  this  way,  regional  delay  in  excitation,  ischemia,  injury,  and 
infarction  could  be  simulated.  This  model  could  be  further  refined  to  model 
subendocardia  1 ,  intramural,  subepicardial ,  or  transmural  ischemia,  injury,  or 
infarction  in  each  of  these  regions.  The  volume  conductor  was  assumed  to  be 
an  infinite  volume  conductor  of  uniform  resistivity.  Recording  sites  in  the 
medium  were  placed  to  represent  the  12-lead  ECG  locations  on  an  adult  male 
torso.  While  the  volume  conductor  was  unrealistic  in  that  internal  torso 
inhomogeneities,  such  as  intracardiac  blood  mass  and  lungs,  and  an  external 
irregular  torso  boundary  with  air  outside  was  not  simulated,  this  was  the 
first  simulation  of  3D  activation  and  recovery  from  a  3D  model  of  the  actual 
heart  geometry.  This  simulation  produced  a  12-lead  ECG  that  mimicked 
regional  ischemia,  injury,  and  infarction  in  a  way  entirely  consistent  with 
clinical  ECGs  recorded  in  these  situations.  ECGs  produced  by  the  simulation 
and  read  by  experienced  el ectrocardiographers  were  interpreted  as 
representing  the  process  simulated. 
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In  1978,  Horan  and  associates  further  refined  the  heart  geometry  by 
dividing  the  ventricles  of  a  real  human  heart  into  3.2  mm3  regions  (143). 
The  volume  conductor  was  again  assumed  to  be  an  infinite  homogeneous  resis¬ 
tive  medium.  They  also  included  activation  times  and  a  base  to  apex  and 
epicardial  to  endocardial  gradient  in  action  potential  durations  that  matched 
published  data.  They  used  a  simplified,  trapezoidal  shape  of  action  poten¬ 
tial  waveshape.  For  computational  simplification,  the  small  regions  were 
combined  together  in  groups  of  8  by  first  calculating  the  differences  between 
the  8.  The  contributions  of  each  group  to  the  extracellular  current  field 
were  estimated  from  the  potential  differences  between  the  regions.  Using 
dipole  and  quadripole  contributions  of  the  local  regions,  this  group  investi¬ 
gated  the  theoretical  validity  of  the  concept  of  "primary"  and  "secondary"  T 
waves.  They  calculated  an  "intrinsic"  T  wave  by  simulating  simultaneous 
excitations  of  all  regions  and  compared  this  generator  output  to  that  calcu¬ 
lated  by  the  method  of  Abildskov  et  al.  They  found  a  gocd  correlation,  but 
not  complete  agreement  between  the  two. 

Miller  and  Geselowitz  added  a  realistic  external  torso  boundary  simula¬ 
tion  and  more  realistic  action  potential  waveshapes  to  compute  body  surface 
ECG  maps  of  the  excitation-recovery  process  (52,53).  They  also  produced  12- 
lead  ECGs  of  simulated  normals,  as  well  as  multiple  combinations  of  regional 
subendocardial,  subepicardial,  and  transmural  ischemia,  injury,  and  infarc¬ 
tion.  This  is  the  most  definitive  simulation  to  date  that  incorporates  meas¬ 
ured  heart  geometry,  activation  sequences,  and  distributed  gradients  of 
recovery  into  a  realistic  bounded  male  torso.  Our  proposed  simulation  will 
refine  the  cardiac  generator  simulation  even  further  (to  a  1  mm 3  grid)  and 
will  include  internal  inhomogeneities,  i.e.,  intracardiac  blood  mass  and 
lungs  into  the  model,  and  will  develop  a  library  of  torso  shapes  and  internal 
torso  geometries  and  properties  such  that  when  completed  the  simulations  can 
be  generalized  to  the  pilot  population  at  large,  both  at  rest  and  with 
exercise.  It  will  also  include  systolic  wall  motion  into  the  simulation. 
When  this  is  accomplished,  we  expect  to  have  incorporated  all  the  first-order 
variables  into  the  simulation  for  both  excitation  and  recovery  (except  ven¬ 
tricular  ectopic  beats  and  bundle  branch  blocks)  that  are  known  to  effect  the 
potential  distributions  on  the  body  surface  at  rest  and  during  exercise. 

In  its  present  form,  the  forward  simulation  includes  the  simulation  of 
fascicular  blocks  in  either  or  both  the  right  or  left  bundle  systems.  This 
allows  for  the  simulation  of  ischemia,  injury,  and  infarction  in  the  presence 
of  left  superior,  middle,  or  inferior  fascicular  block.  This  is  accomplished 
in  the  model  by  eliminating  start  points  in  each  of  these  distributions  where 
the  left  fascicular  system  inserts  into  the  Purkinje  network.  At  some  point 
in  the  development  of  this  forward  model,  it  will  be  important  to  include  the 
development  of  a  simulation  of  complete  left  and  right  bundle  branch  block 
into  the  system.  The  ability  to  predict  the  presence  of  ischemia,  injury  and 
infarction,  and/or  hypertrophy  in  the  presence  of  these  major  conduction 
abnormal i ties  is  of  considerable  interest,  and  the  timing  of  its  development 
will  be  dictated  in  this  effort  by  the  relative  significance  of  this  subset 
of  the  problem  in  the  Air  Force  pilot  population.  Specifically  at  some  point 
if  the  results  from  forward  simulation  effort  have  led  us  to  the  place  where 
the  sensitivity  and  specificity  for  the  detection  of  asymptomatic  coronary 
disease  is  such  that  the  case  finding  approaches  95%  in  subjects  without 
bundle  branch  block,  then  it  would  seem  prudent  to  include  these  conduction 
defects  in  the  model.  At  that  time  this  group  of  subjects  would  begin  to 


represent  a  significant  proportion  of  the  residual  pilot  population  that 
might  have  undetected  asymptomatic  coronary  artery  disease.  This  will  be 
discussed  in  more  detail  along  with  the  inclusion  of  regional  intramural 
slowing  of  conduction  and  anisotropies  as  relates  to  ischemia  and  infarction 
in  the  section  on  "Numerical  Experiments  with  the  Forward  Model  and  Sequence 
of  Its  Development". 


Torso  Transfer  Impedance  Changes  at  Rest  and  Exercise 


Ventricular  Wall  Motion  with  Systole 


In  developing  the  simulation  of  depolarization,  the  motion  of  the  heart 
in  man  could  be  ignored  since  ventricular  excitation  occurs  during  the  quiet 
end  diastolic  phase  (or  immediately  presystole).  In  the  modelling  of  the 
recovery  phase  or  ST-T  segment  of  the  ECG,  which  occurs  throughout  systole, 
the  motion  of  the  heart  during  this  time  must  be  considered  in  the  formula¬ 
tion  of  the  transfer  matrices  between  local  heart  regions  and  the  body  sur¬ 
face.  Since  potential  differences  between  2  points  in  an  infinite  homoge¬ 
neous  resistive  volume  conductor  vary  as  the  square  of  the  distance  between 
these  2  points,  it  follows  that  variations  in  potential  at  the  surface  of  the 
thorax  of  man  generated  from  cardiac  dipole  current  generators  are  quite 
sensitive  to  changes  in  position  of  the  generators  with  systole,  especially 
in  those  generators  close  to  the  torso  surface.  Potential  changes  on  the  body 
surface  produced  from  those  regions  deep  in  the  chest  will  likewise  be  less 
sensitive  to  systolic  wall  motion. 

Cardiac  size  and  output  measurements  in  man  during  rest  and  exercise  have 
been  studied  using  simple  hemodynamic  and  radiographic  means  by  Meek,  McCrea, 
and  Eyster  (144,145),  and  by  others  including  early  work  by  Rushmer's  group 
and  Sarnoff  and  Mitchell  (146-151).  These  data  included  not  only  the  con¬ 
traction  and  volume  changes  of  the  ventricles  at  rest  and  with  exercise,  but 
translational  effects.  Conflicting  data  with  regard  to  a  "wringing"  type  of 
rotational  motion  around  the  long  axis  was  described  in  these  early  papers. 
The  base  of  the  left  ventricle  was  generally  seen  to  rotate  clockwise  around 
the  long  axis,  and  the  apex  rotated  counterclockwise.  The  hemodynamic  and 
volumetric  effects  of  respiration  were  also  delineated  in  these  early 
reports. 


These  observati 
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ons  have  been  refined  in  a  semiquantitative  way  in  dogs  by 
(146),  and  by  Rushmer  and  collaborators  (148)  by  suturing 
the  surface  of  the  ventricles  and  observing  them  fluoro- 
al lowed  for  more  detailed  assessment  of  regional  wall 
ics,  beat  by  beat  in  the  conscious  animal  over  extended 
with  various  interventions.  Beginning  in  1963,  Harrison, 
Id,  Glick,  and  Mason  (152-154)  in  Bethesda,  McDonald  (155) 
ngles,  Daughters,  Stinson,  and  Alderman  (156)  in  Palo  Alto 
studies  in  man  with  the  use  of  small  epicardial  or  intra- 
rs  implanted  at  the  time  of  open  heart  surgery.  Cardiac 
in  the  intact  conscious  man  could  now  be  done  in  compari- 
1  biplane  contrast  ventriculograms,  and  once  validated, 
cineradiographic  system,  serial  measurements  can  be  made 


to  determine  the  effects  of  various  physiological  and  pharmacological 
interventions  on  epicardial  dynamics. 

Wildenthal  in  1969  (157)  demonstrated  in  experimental  animals  significant 
differences  in  beat-to-beat  ventricular  dynamics  measured  from  epicardial  and 
midwall  markers.  The  Ingels  paper  (156)  with  tantalum  markers  in  the  midwall 
is  the  most  relevant  to  the  simulation  of  the  ECG  during  the  recovery  phase, 
since  the  present  version  of  the  computer  simulation  coalesces  all  current 
generators  in  a  given  segment  to  a  dipole  at  the  centroid  of  each  of  the 
myocardial  segments.  Figure  10  shows  the  local  rotational  and  translational 
effects  (as  demonstrated  especially  by  Ingles  et  al.)  in  each  of  the  12  left 
ventricular  segments  that  will  be  used  in  the  simulation  of  systole  at  rest. 


apical 


Figure  10.  Examples  of  the  complex  positional  changes  from  end  diastole  to 
peak  systole  are  shown  here  for  a  few  illustrative  segments.  The 
basal  segments  will  move  down,  inward,  and  rotate  clockwise  with 
reference  to  the  long  axis  as  viewed  from  the  apex.  Apical  seg¬ 
ments  will  move  upward,  inward,  and  rotate  counterclockwise.  Mid 
wall  segments  will  move  mainly  inward. 


Coronary  artery  bifurcations  were  used  by  Kong,  Morris  and  McIntosh  (158) 
to  assess  regional  wall  function  in  patients  with  coronary  artery  disease. 
This  method  was  checked  against  the  method  using  epicardial  metal  markers  in 
the  experimental  laboratory  and  it  was  found  that  for  the  first  few  beats 
after  a  coronary  injection  regional  wall  motion  of  the  epicardial  markers  was 
not  affected.  This  method  was  labor  intensive  but  did  demonstrate  in  this 


small  series  that  there  was  systolic  lengthening  on  the  distal  distribution 
of  high  grade  stenotic  vessels,  and  decreased  and  delayed  or  absent  regional 
shortening  in  infarcted  regions  distal  to  complete  coronary  obstructions.  It 
will  be  important  to  incorporate  these  observations  into  the  simulation  of 
regional  ischemia,  and  infarction. 

The  extensive  recent  literature  in  regional  wall  motion  changes  by 
nuclear  angiographic  methods  and  2D  Echo  is  pertinent  in  this  regard.  The 
earliest  marker  of  regional  ischemia  has  been  shown  to  be  thinning  of  the 
left  ventricular  wall,  and  slight  outward  bulging.  Transmural  infarction  is 
associated  with  a  thick  collagen  scar,  and  local  akinesis.  These  effects, 
and  the  well -documented  translocational  effects  of  inferior  infarction  from 
right  coronary  occlusion,  in  particular,  will  need  to  be  incorporated  into 
the  location  changes  and  transfer  matrix  changes.  Criteria  for  which  of 
these  variables  need  to  be  included  will  result  from  the  numerical  analysis 
studies  and  the  numerical  experiments  referred  to  in  the  section  "Numerical 
Experiments  with  the  Model  and  Sequence  of  Its  Development"  (pp.  69-72). 


Exercise-Induced  Changes 

A  reasonable  body  of  older  literature  deals  with  the  changes  in  the 
cardiac  size  and  output  with  exercise  (144-149).  These  studies  used  standard 
biplane  radiographs,  and  included  atria  and  great  vessels  in  general  in  the 
measurement  of  heart  volumes  and  dimension  changes  with  exercise.  The  mea¬ 
surement  of  the  left  ventricular  wall  regional  changes  and  motion  with  near 
maximum  effort  is  mentioned  by  Ingels  and  coworkers  (156)  as  one  of  the 
potential  uses  of  the  midwall  implanted  marker  technique,  but  no  such  speci¬ 
fic  data  was  reported.  Very  limited  data  was  reported  by  Braunwald, 
Goldblatt,  Harrison,  and  Mason  (154),  looking  at  left  ventricular  length  by 
means  of  2  epicardial  markers  (one  near  the  apex,  the  other  near  the  base)  in 
five  subjects  during  mild  exercise.  They  found  little  change  (average  6.5% 
decrease)  in  this  measurement  with  light  exercise.  They  concluded  from 
review  of  the  older  data  mentioned  above  and  their  own  measurements  that: 
1)  there  was  a  modest  decrease  in  left  ventricular  end  diastolic  and  end 
systolic  volume  with  exercise;  2)  the  main  change  with  exercise  was  a  signi¬ 
ficant  increase  in  contractility,  most  likely  due  to  increased  sympathetic 
tone;  and  3)  the  change  in  cardiac  output  was  the  result  of  an  increase  in 
systolic  emptying  and  an  increase  in  heart  rate. 

The  data  needed  to  simulate  right  and  left  ventricular  dimensional 
changes  with  exercise  in  presumably  normal  adult  men  is  likely  present  in  the 
raw  data  of  some  of  the  2D  Echo  data  now  being  reported,  but  the  published 
reports  are  concerned  with  regional  and  global  changes  in  dynamics  and  are 
generally  corrected  for  the  movement  of  the  heart  downward  by  normalizing  the 
systolic  frames  used  to  the  aortic  root.  Most  reports  also  use  the  subcostal 
long  axis  view  where  the  distance  to  the  overlying  thorax  is  not  measured. 
Nevertheless,  in  reviewing  this  problem  with  Ginston  (159),  who  has  reported 
on  the  2D  Echo  changes  with  exercise  (160),  it  is  evident  that  in  echogenic 
subjects,  the  nature  of  these  local  dimensional  changes  can  be  assessed,  and 
incorporated  into  the  model  if  the  early  numerical  experiments  and  resolution 
studies  indicate  that  this  is  needed.  This  will  require  fixing  the 
relationship  between  the  thorax  of  the  subject  and  the  transducer,  either  by 


using  the  apical  long  axis  views,  or  parasternal  views  at  rest  and  during 
exercise. 

Regarding  the  changes  in  torso  resistivities  with  exercise,  detailed  data 
of  the  kind  needed  in  the  simulation  is  not  available  in  the  current  litera¬ 
ture,  and  will  need  to  be  generated  as  a  part  of  the  ongoing  effort  in  later 
years  of  this  project.  The  measurement  of  torso  impedance  changes  has  been 
done  as  a  part  of  an  earlier  study  reported  by  us  (161)  in  the  baboon,  and 
preliminary  studies  in  volunteer  subjects  indicate  the  feasibility  of  making 
these  measurements  in  man.  Transthoracic  impedance  measurements,  as  reported 
by  Krohn  (162),  are  done  routinely  as  a  qualitative  measure  of  cardiac  output 
trends,  using  commercially  available  equipment.  The  objective  of  the  work 
being  proposed  by  us  will  be  to  make  these  measurements  quantitative  by 
measuring  small  currents  injected  at  one  electrode  site  of  a  limited  map 
array  at  all  other  sites,  and  repeating  this  at  each  electrode  until  an 
impedance  map  is  generated. 


Optimal  Electrode  Number  and  Location 

Since  the  early  work  by  Taccardi  (163-166),  and  by  Horan,  Flowers,  and 
Brody  (167),  on  ECG  body  surface  mapping  of  heart  potential  distribution, 
there  has  been  steady  progress  in  establishing  regional  ECG  surface  map 
changes  that  relate  to  regional  infarct  changes  (168-170).  Horan,  Flowers, 
and  Brody  reported  early  work  (171-173)  on  principal  factor  waveforms  extrac¬ 
tion  and  multipolar  content  of  the  map  data.  Using  similar  numerical 
(empiric)  methods,  Horan  and  coworkers  (174),  Lux,  Burgess  and  coworkers 
(175,176),  Barr,  Spach  and  coworkers  (177-179),  and  Kornreich  and  coworkers 
(180,181),  have  shown  that  the  total  ECG  surface  map  from  150-200  points  can 
be  reproduced  to  some  acceptable  resolution  with  15  to  32  leads,  depending  on 
the  diagnostic  information  content  and  the  resolution  desired.  In  some 
instances  the  resolution  needed  for  the  present  effort  was  not  attempted. 


Using  the  total  body  surface  map  simulation,  which  with  our  adaptation  of 
the  Gerlernter  Swihart  inhomogeneous  torso  simulation  generated  a  1874-point 
ECG  surface  map,  we  looked  at  the  map  potential  distribution  from  unit 
dipoles  at  each  of  the  local  left  ventricular  segments  (182).  Predictably 
some  segments,  like  the  anteroapical  ones,  produced  large  peak  potentials 
with  well  localized  distributions  on  the  body  surface,  and  others  like  the 


posterobasal  segment  produced  a  much  smaller  peak  surface  distributed  over  a 
much  wider  area.  Generally  speaking,  however,  there  were  unique  surface 
distributions  for  each  of  the  major  left  ventricular  segments.  These 
topographical  relationships  for  unit  dipoles  are  shown  in  Figures  11  and  12. 


Using  this  anatomical  and  physiologically  oriented  information,  we  plan 
to  use  the  simulation  to  define  the  expected  resolution  for  regional  ischemia 
in  each  of  the  12  segments  of  our  current  version  of  the  simulation.  Speci¬ 
fically  how  much  ischemia  does  it  take  in  each  of  the  local  segments  to 
produce  surface  ECG  changes  above  the  noise  level  expected  in  exercise  ECGs? 
Where  does  this  appear  on  the  body  surface?  What  specific  criteria  will  be 
needed  to  detect  this  degree  of  ischemia?  From  the  distribution  of  unit 
dipoles  in  Figure  11  it  is  clear  that  similar  degrees  of  subendocardial 
ischemia  in  the  apical  segments  will  be  detected  with,  for  example,  200-300 
V  of  ST  segment  change  in  anterior  precordial  leads  whereas  that  degree  of 
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Figure  11.  The  flat  surface  represents  the  torso  surface  opened  vertically  at 
the  right  mid  axillary  line  with  the  back  unrolled  to  the  left  on 
the  diagram.  Each  left  ventricular  unit  dipole  is  clipped  at  50% 
of  positive  potential  to  show  the  regional  localization  of  that 
segment  on  the  body  surface  and  the  strength  of  each  local  unit 
dipole  on  the  surface. 
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Figure  12.  In  this  illustration  the  surface  distribution  of  all  the  left 
ventricular  unit  dipoles  shown  individually  in  Figure  11  are 
projected  onto  the  surface.  The  coronary  artery  supply  to  each 
segment  is  shown  at  the  top.  It  will  be  noted  in  the  bottom  frame 
that  the  local  distribution  of  each  of  the  segments  is  well  pre¬ 
served.  The  height  of  each  segment  distribution  shown  is  related 
to  the  sensitivity  of  the  ECG  surface  map  to  the  local  left 
ventricular  region. 


ischemia  in  the  posterobasal  segment  would  produce  only  about  20-30  /jV  of  ST 
change  localized  preferentially  to  leads  on  the  back.  With  existing  signal 
processing  techniques,  this  degree  of  change  should  be  well  within  the  limits 
of  the  technology. 

The  strength  of  this  approach  is  that  the  model  is  put  together  from  a 
1  mm3  "cell"  to  the  body  surface  from  anatomical  and  physiological  informa¬ 
tion.  If  the  simulation  contains  all  the  first-order  variables  that  influ¬ 
ence  the  output,  then  changes  seen  on  the  surface  have  a  direct  anatomical 
and  physiological  base.  The  changes  so  observed  yield  directly  to  anatomical 
and  physiological  diagnostic  statements.  On  the  other  hand,  empiric  methods 
such  as  multipolar  content  and  multipolar  reductions  of  surface  data,  while 
computat ional ly  more  efficient  and  mathematical ly  more  direct,  do  not  have  a 
direct  analogy  to  the  anatomy  of  the  heart.  For  example,  some  quadrapole  or 
octapole  term  may  be  quite  outside  two  standard  deviations  of  that  seen  in  a 
significant  population  of  normal  subjects,  but  large  clinical  correlative 
studies  will  need  to  be  done  to  establish  the  physiological  or  anatomical 
significance  of  the  specific  coefficients  that  were  found  to  be  abnormal. 

It  also  follows  that  maps  generated  by  the  Lux  leads,  and  Eigenvector 
algorithms,  or  by  the  Barr-Spach  24-lead  set  can  be  produced  from  the  larger 
lead  set  maps  that  are  generated  by  various  combinations  of  regional  ischemia 
and  infarction  using  our  simulation.  The  ability  of  each  of  these  data 
reduction  techniques  to  reproduce  the  surface  ECG  map  to  the  resolution 
needed  to  detect  the  specific  change  can  be  directly  assessed.  The  simula¬ 
tion  can  be  used  as  a  tool  to  develop  a  theoretical  base  for  the  interpreta¬ 
tion  of  various  abnormal i ties  seen  in  the  limited  set,  and  to  explore  whether 
diagnostically  important  areas  on  the  body  surface  need  lower,  or  higher 
concentration  of  leads  to  achieve  the  needed  sensitivity  and  specificity  for 
either  ischemia  or  infarction. 
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FORWARD  MODEL  DESIGN  AND  DEVELOPMENT 


Task  1  has  been  completed  as  of  December  31,  1984.  The  review  of  the 
literature,  and  the  design  and  upgrading  of  the  Sel vester-Sol omon  fine  grid 
(1  mm3)  forward  propagation  model  of  human  ventricular  excitation  to  include 
repolarization  have  been  accomplished.  The  presentation  of  this  document,  a 
USAFSAM  technical  report,  summarizes  these  findings  and  completes  Task  1. 
The  design  of  the  upgraded  model  is  presented  as  a  Treatise  beginning  in  the 
Appendix.  Its  development  and  implementation  are  described  as  Task  2. 


Task  2:  The  Forward  Model  Development  has  two  major  subtasks  (2A  and  2B) 
that,  because  of  limited  funding,  will  be  run  serially  (2-year  effort). 
Subtask  2A  (a  1-year  effort)  is  the  actual  development  of  a  working  version 
of  the  model  with  one  or  more  torso  simulations  representative  of  the  Air 
Force  pilot  population  and  typical  examples  of  local  ischemia,  injury,  and 
infarct.  Significant  programmatic  efficiencies  described  next  would  result 
if  the  funding  can  be  upgraded  at  the  end  of  2A  to  allow  Subtask  2B  and  Task 
3  to  develop  concurrently.  Subtask  2B  (a  1-year  effort  at  the  current  fund¬ 
ing  level)  is  the  development  of  the  table-driven  subroutines  to  modify  the 
geometry,  activation,  and  action  potentials  in  order  to  simulate  at  will  a 
comprehensive  catalog  of  regional  ischemia,  injury  and  infarction,  at  rest 
and  with  exercise.  Upon  the  completion  of  Task  2  (in  two  years  at  current 
funding  levels)  we  expect  to  begin  Task  3  (Optimal  Leads  and  Criteria)  along 
with  early  work  on  the  validation  studies. 

Subtask  2A  is  the  actual  development  of  a  working  version  of  the  model 
with  a  typical  torso  and  a  few  typical  areas  of  ischemia,  injury,  and  infarct 
(see  "Flow  Chart  of  the  Forward  Model  Program").  In  addition  to  the  Forward 
Computer  Model  development  described  in  detail,  this  includes: 

a.  Torso  tank  verification  of  the  Gelernter-Swihart  simulation  for 
a  broad  range,  representative  set  of  local  dipoles  and  torso  resistiv¬ 
ities  to  include  the  expected  extrema  of  locations,  torso  shapes,  and 
resistivities. 

b.  Torso  tank  studies  on  the  potential  errors  in  lead  placement, 
i nterel ectrode  res i sti vi ty  measurements ,  and  methods  for  verifying  the 
relationship  between  the  leads. 

c.  Measurement  of  torso  geometry  and  impedance  changes  in  man  with 
near  maximal  and  maximal  exercise. 

This  effort  requires  1  man  year  of  physicist/analyst  time  and  1/4  man  year  of 
cardiologist-physiologist  time. 

Subtask  2B  is  a  separate  effort  to  develop  the  table-driven  subroutines 
that  modify  the  local  initial  conditions  of  size,  transmural  extent,  and 
location  of  ischemia,  injury,  and  infarct,  and  whether  these  are  single  or 
multiple.  These  modifying  subroutines  will  be  able  to  change  heart  location 
and  long  axis  orientation  in  the  chest,  etc.,  simulate  local  or  generalized 
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hypertrophy  or  dilatation  of  either  or  both  ventricles,  and  simulate  varia¬ 
tions  in  body  build  and  resistivities.  Torso  properties,  heart  locations, 
etc.,  may  be  modified  from  on  line  torso  resistivity,  anthropometric  (skin- 
folds,  torso  geometry,  etc.),  and  2-D  echocardiography  measurements,  if  error 
analysis  and  analysis  of  the  causes  of  variances  indicate  this  may  improve 
the  resolution  of  the  program.  This  subtask  requires  1  man  year  of  a  physi¬ 
cist/analyst/programmer  and  1/4  man  year  of  cardiologist/physiologist  time. 

By  running  Subtasks  28  and  Task  3  concurrently,  significant  efficiencies 
and  facilitation  of  the  overall  objectives  will  be  realized.  The  following 
staffing  will  be  required:  1  man  year  of  a  physicist/analyst,  1  man  year  of 
an  experienced  scientific  utility  programmer,  1  man  year  of  a  research 
technician,  and  1/2  man  year  of  cardiologist/physiologist  time. 

IN  SUMMARY:  Our  forward  model  design  begins  with  the  simulation  of  normal 
heart  depolarization  and  repolarization  in  a  typical  inhomogeneous  male 
torso.  This  model  is  an  extension  of  the  work  we  have  already  done  on 
simulating  the  QRS  activation  and  will  include  the  physiology  of  action 
potentials,  particularly  the  epicardial  to  endocardial  and  base  to  apex 
gradients  of  action  potential  duration  and  downslope,  and  the  geometry  of  the 
heart  at  a  1  mm3  grid.  This  normal  simulation  will  be  modified  to  include 
the  geometry  of  infarct  and  ischemia  as  well  as  the  motion  of  the  walls 
(dipoles)  with  contraction.  At  this  point,  we  will  be  able  to  test  two 
hypotheses:  1)  the  Mil ler-Geselowitz  model  of  ischemia  where  there  is  a 
series  of  border  zones  surrounding  the  dead  tissue  with  progressively  less 
ischemia  the  farther  from  the  necrotic  area  one  measures;  and  2)  that  there 
is  a  sharp  border  separating  infarcted  tissue  from  healthy  tissue  and  it  is 
the  complexity  of  the  geometry  of  the  infarct  that  produces  the  observed 
electrocardiographic  changes. 

Once  this  is  established,  we  will  use  a  1  mm3  descretization  of  the  heart 
as  a  starting  point.  Since  developing  transfer  impedances  from  such  a  large 
number  of  points  is  unwieldy,  we  will  group  these  points  into  12  subsets  for 
the  left  ventricle  and  8  for  the  right  ventricle  for  a  total  of  20  dipoles. 
We  will  check  the  Thiry-Rosenberg  method  of  representing  the  cardiac  dipoles 
against  this  1-mm  grid  model  with  20  dipoles  to  see  which  method  is  opera¬ 
tionally  most  tractable  and  computationally  the  most  economical.  After  the 
optimal  subsets  have  been  determined,  we  will  develop  transfer  impedances  for 
a  variety  of  normal  male  torsos.  The  next  step  is  to  generate  ECG  maps.  We 
will  be  able  to  display  scalars,  12-lead,  VCG's,  perspective  plots,  isopoten¬ 
tial  plots  (in  color  or  black  and  white)  with  various  sized  single  or  mul¬ 
tiple  areas  of  infarct,  injury,  and  ischemia  in  normal  hearts  or  associated 
with  hypertrophy,  dilatation,  and/or  conduction  abnormalities.  Isoarea  QRS, 
ST,  T,  and  QRST  maps  will  also  be  examined. 
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Detailed  Description  of  the  Forward 
Computer  Model  Development 


The  Task  2  model  development  consists  of  writing  the  source  code,  debug¬ 
ging,  testing,  and  documenting  the  three  programs  that  constitute  the  forward 
simulation. 

1.  The  wave  propagation  program  simulates  the  myocardial  activation 
front. 

2.  The  equivalent  cardiac  source  program  generates  equivalent 
myocardial  sources  for  depolarization  repolarization. 

3.  The  ECG  program  computes  body  surface  ECG  data  from  the  equivalent 
cardiac  generators  and  the  inhomogeneous  torso  transfer  impedances. 

Subtask  2A-1:  Propagation  Model. 

We  will  do  the  following: 

1.  Create  a  normal  human  heart  and  Purkinje  geometry  file  for  input. 

2.  Write  the  propagation  program  source  code  to  execute  on  the  VAX- 
7 80/VMS  system 

3.  Generate  an  output  file  of  the  myocardial  activation  data  for  the 
normal  human  heart  that  is  format  compatible  with  the  equivalent  cardiac 
source  program. 

4.  Run  the  propagation  program  on  a  spherical  geometry  and  verify 
against  the  analytical  activation. 

5.  Produce  graphic  output  in  CRT  and  hard  copy  of  the  activation  in  the 
normal  heart  geometry  at  various  levels  of  vertical  and  horizontal  cross 
section  to  verify  that  the  programs'  algorithms  are  error  free. 

6.  Document  the  source  code  and  describe  the  input/output  data  and  file 
formats. 

Subtask  2A-2:  Equivalent  Source  Program. 

This  program  includes  the  1  mm^  descretization  and  the  grouping  of  these 
into  12  LV  and  8  RV  myocardial  segments  (20  dipoles).  We  will  do  the 
following: 

1.  Create  a  nominal  value  file  for  the  action  potential  profiles 
throughout  the  myocardium. 

2.  Write  the  source  code  for  the  equivalent  source  program  to  run  on 
the  VAX-780/VMS  system. 

3.  Verify  that  the  program  will  read  the  compatible  activation  data 
from  the  wave  propagation  program  under  all  possible  variations. 
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4.  Have  the  program  generate  an  output  file  of  the  equivalent  cardiac 
source  that  is  compatible  with  the  ECG  program. 


5.  Display  graphically  the  equivalent  cardiac  generators  (action  poten¬ 
tial  profiles)  on  CRT  and  hard  copy  and  verify  the  program  output  by 
comparison  with  hand  computed  values  based  on  the  input  data  files. 

6.  Test  the  program  with  a  10  mm  cube  of  simulated  myocardium  and 
nominal  action  potential  profiles. 

7.  Display  graphically  (CRT  and  hard  copy)  the  output  of  the  20  dipole 
myocardial  segment  subdivision. 

8.  Document  the  source  code  and  describe  the  input/output  data  and 
formats. 

Subtask  2A-3:  ECG  Data  Program. 

We  will  do  the  following: 

1.  Solve  the  inhomogeneous  boundary  value  problem  for  a  normal  set  of 
transfer  impedances  from  heart  generators  to  torso  surface  leads. 

2.  Write  source  code  to  compute  ECG  data  from  the  input  files. 

3.  Verify  that  the  source  code  can  read  the  equivalent  cardiac  source 
input  file  under  all  possible  variations. 

4.  Code  the  ECG  data  program  to  generate  an  output  file  formatted  in 
sections  which  contain  contiguous  blocks  of  data  for  each  type  of  output, 
such  as  12  lead  ECG,  McFee  vector,  Frank  vector,  body  surface  map,  etc. 

5.  Test  the  boundary  value  solution  against  the  electrolytic  torso  tank 
for  several  sources  to  demonstrate  the  agreement  between  the  mathematical 
solution  and  physical  measurements. 

6.  Graphically  display  the  ECG  output  data  file  for  unit  transfer 
impedances  and  unit  orthogonal  dipole  sources  to  verify  the  program 
al gorithms. 

7.  Document  the  source  code  and  describe  the  input/output  data  and 
formats. 


Studies  with  the  Forward  Computer  Model 

Subtask  2A-4.  Verification  of  the  Transfer  Impedances  on  the  Electrolytic 
Tank: 

Since  the  transfer  impedances  of  the  forward  model  simulation  are  derived 
from  a  mathematical  solution  to  a  boundary  value  problem  using  matrix  and 
iterative  methods,  it  becomes  necessary  to  verify  the  mathematical  solutions 
in  an  experimental  configuration.  Verification  will  consist  of  computing  the 
transfer  impedances  for  the  electrolytic  torso  tank  that  contains  a  number  of 


dipole  sources  and  comparing  them  with  the  actual  measured  values.  Results 
of  such  a  verification  experiment  determine  the  errors  in  the  boundary  value 
solution  relative  to  the  location  and  direction  of  the  dipoles  imbedded  in 
the  volume  conductor.  There  is  no  other  reliable  test  of  the  program  code, 
convergence,  and  accuracy  than  a  physical  verification. 

Subtask  2A-5.  Electrode  Misplacement  and  Reconstruction: 

An  experimental  study  will  be  conducted  on  the  electrolytic  tank  to 
develop  a  method  of  identifying  chest  lead  misplacements  and  reconstruct  the 
correct  ECG.  This  procedure  will  investigate  the  use  of  small  sinusoidal 
currents  impressed  on  the  torso  to  detect  and  correct  for  improper  transfer 
impedances.  Safety  considerations  for  use  on  humans  will  be  identified  and 
establ ished. 

Subtask  2A-6.  Effects  of  Exercise  on  Torso  Transfer  Impedances: 

A  major  part  of  the  forward  model  application  is  in  determining  optimal 
leads  and  criteria  at  exercise.  For  this  reason,  the  torso  volume  conductor 
changes  that  occur  in  normal  adult  males  at  exercise  will  be  measured. 
Changes  in  volume  conductor  characteristics  measured  between  rest  and  exer¬ 
cise  will  be  computed  with  the  boundary  value  solution  program  to  estimate 
the  properties  of  the  torso  transfer  impedances  at  exercise.  These  effects 
will  become  a  part  of  the  variations  included  in  the  transfer  impedance 
library  for  the  forward  simulation  program. 

Subtask  2B:  Develop  the  Input/Output  Data  Generation  Programs: 

We  will  develop  the  input/output  programs  that  produce  the  disk  files 
used  for  input  to  the  forward  ECG  simulation  programs.  These  data  generation 
programs  interact  with  the  operator  and  the  simulation  programs.  An  operator 
responds  to  this  set  of  programs  to  generate  various  input  conditions  that 
produce  a  desired  ECG  simulation.  The  four  programs  to  be  developed  to 
interface  the  operator  with  the  simulation  programs  are  as  follows: 

1.  Heart  geometry  alteration  routine. 

2.  Purkinje  modification  including  starting  points. 

3.  Action  potential  distribution  program  for  specifying  the  regional 

distribution  parameters  for  equivalent  cardiac  sources. 

4.  Transfer  impedance  selection  program  that  retrieves  a  specific  set 

of  transfer  impedances  from  a  cataloged  library. 

Task  2C:  Torso  Transfer  Impedance  Library: 

The  Gelernter-Swihart  method  of  computing  a  matrix  of  torso  impedances 
requires  several  hours  of  time  for  each  variation  of  heart  location  or  geome¬ 
try,  torso  shape,  or  resistivity.  The  development  of  a  comprehensive  library 
of  torso  transfer  impedances  will  be  a  major  ongoing  effort  requiring  1-3 
years  for  completion  depending  on  the  availability  of  computer  time  and 
funds.  At  the  start  up  of  Task  2,  we  will  immediately  begin  generation  of 
the  library  of  inhomogeneous  torso  (volume  conductor)  transfer  impedances 


from  equivalent  source  locations  in  the  heart  to  torso  surface  electrode 
sites. 

Initially  a  20-dipole  equivalent  heart  source  will  be  used  to  represent 
the  electrical  activity,  and  the  transfer  impedances  will  be  determined  for 
the  three  orthogonal  components  of  each  of  the  20  dipoles.  Transfer  impe¬ 
dances  will  be  computed  for  numerous  variations  about  the  nominal  (typical) 
case.  These  variations  include  short,  medium,  and  tall  torso  dimensions 
(relative  to  Air  Force  standards),  as  well  as  thin,  typical,  and  heavy  body 
weights. 

Anterior-posterior  to  lateral  ratios  will  be  included.  Another  first- 
order  transfer  impedance  variable  is  the  position  and  orientation  of  the 
heart  in  the  chest.  Transfer  impedances  will  be  computed  for  various  aortic 
root  locations  and  Euler  angles  specifying  the  orientation  of  the  heart. 
Initially,  three  sets  of  transfer  impedances  for  each  distinct  torso  condi¬ 
tion  will  be  computed;  one  each  in  diastole,  mid  systole,  and  end  systole. 
Interpolation  through  these  three  cases  will  be  utilized  to  obtain  transfer 
impedances  for  ST-T  wave  computation  while  the  heart  is  in  motion. 

Variations  in  lung  and  blood  mass  conductivities  in  the  Air  Force  pilot 
population  is  likely  to  be  a  second-order  effect  and  may  not  need  to  be 
included.  However,  such  variations  would  be  included  if  more  detailed 
analysis  of  the  subjects  indicates  that  it  is  required. 


Val idation 

Autopsy:  The  standard  12-lead  ECG  taken  no  more  than  10  days  before 
death  will  be  simulated  in  the  following  6  male  patients,  age  25  through  55. 
Two  subjects  will  have  no  organic  heart  disease  on  post  mortem  examination  by 
the  quantitative  planometric  method  of  Idecker  and  Hackel:  two  cases  will 
have  at  least  one  coronary  artery  showing  at  least  50%  cross-sectional 
luminal  narrowing  and  no  evidence  of  myocardial  infarct,  and  two  subjects 
will  have  had  at  least  one  old  infarct.  No  subject  will  have  had  a  recent 
infarct,  myocarditis,  cardiomyopathy,  congenital  or  valvular  heart  disease, 
evidence  for  passive  congestion  or  abnormal  fluid  accumulation  in  any  body 
cavity,  or  any  electrolyte  abnormality  at  the  time  the  ECG  was  taken.  Cases 
will  be  selected  that  are  of  the  same  general  height,  weight,  body  build  and 
torso  configurations  as  found  in  typical  airmen.  The  12-lead  ECG  of  all  6 
subjects  will  be  simulated  and  expected  to  pass  the  following  performance 
tests:  when  compared  to  the  last  antemortem  ECG,  the  simulated  ECG  shall 
have  1)  less  than  10  degrees  difference  in  frontal  plane  mean  axis  of  both 
QRS  and  T  wave;  2)  less  than  10-msec  difference  in  both  the  QRS  interval  and 
the  interval  from  the  onset  of  QRS  to  the  peak  of  the  T  wave  in  those  leads 
with  unambiguous  fiducial  points;  and  3)  less  than  100-juv  difference  in  the 
amplitude  of  the  following:  Q,  R,  S,  R',  S',  ST  80  msec  after  J,  and  peak  T. 

Clinical  validation:  The  autopsy  validations  above  establish  the  physi¬ 
ologic  and  anatomic  reasonableness  of  the  simulation.  The  real  utility  of 
the  model,  however,  will  be  in  the  development  of  optimal  lead  sets  and 
criteria  in  those  leads  for  detection  of  ischemia.  The  testing  of  these 
criteria  in  clinical  populations  would  be  the  next  and  manditory  step  in  the 
clinical  validation  of  the  model.  This  is  Task  3  of  this  overall  effort,  and 
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not  dealt  with  directly  in  the  present  document.  In  Task  3,  after  optimal 
leads  and  criteria  have  been  developed,  specific  protocols  will  need  to  be 
set  up  to  test  the  specificity,  sensitivity,  and  predictive  value  of  the 
criteria.  From  a  review  of  the  literature  in  body  surface  ECG  mapping  with 
exercise,  and  from  the  data  summarized  in  this  report,  it  is  reasonable  to 
expect  sensitivities,  and  specificities  in  the  95-98%  range  for  the  detection 
of  unsuspected  ischemia  in  the  asymptomatic  pilot  population.  If  funding  can 
be  increased  in  the  02  yr  to  run  Task  3  concurrently  with  Task  2B,  detailed 
clinical  validation  protocols  that  have  been  previously  developed  will  be 
reviewed  at  that  time. 


Flow  Chart  of  the  Forward  Model  Program 

OVERVIEW:  The  forward  simulation  program  of  human  ECG  data,  as  flow 
charted  here,  is  a  structured  program  executed  from  the  top  down  as  single 
pass,  non-recursive  and  non-reentrant  code.  Control  of  variation  of  input 
conditions  is  through  question  and  answer  dialog.  Simulation  programs  exe¬ 
cute  using  only  input  data  contained  in  current  data  files,  all  ancillary  or 
modification  programs  only  update  these  current  files.  Hence,  input  data  or 
other  operator  inputs  never  alter  the  manner  in  which  any  of  the  simulation 
programs  execute.  Special  programs  are  invoked  by  operator  responses  that 
allow  current  input  data  files  to  be  modified  and  in  some  cases  require 
specific  peripheral  devices  such  as  digitizers  and  plotters.  The  four  current 
input  data  files  required  to  execute  the  entire  simulation  process  are 
discussed  next: 

1.  Current  Heart  Geometry  File:  Normal  cardiac  muscle  and  regions  of 
infarction  and/or  hypertrophy  are  digitized  from  photographs,  drawings,  or 
plots  with  a  digitizing  table,  raster  scanner,  or  similar  device.  This  file 
is  an  input  data  file  for  the  Wave  Propagation  Program  and  is  updated  by  the 
Geometry  Alteration  Program.  A  normal  adult  male  heart  geometry  is  available 
for  normal  ECG  simulations  and  for  testing  purposes.  This  standard  heart 
geometry  can  be  copied  to  the  current  heart  geometry  file  on  program 
initial ization. 

2.  Current  Purkinje  Network  File:  A  standard  set  of  Purkinje  points 
distributed  over  the  endocardium  of  the  standard  heart  geometry,  along  with 
normal  starting  points  are  used  to  initialize  the  current  Purkinje  Network 
File.  The  Purkinje  Modification  Program  is  provided  to  simulate  changes  in 
the  Purkinje  network  and  different  starting  points. 

3.  Current  Action  Potential  Distribution  File:  A  normal  set  of  para¬ 
meters  characterizing  the  action  potential  profiles  in  each  region  of  the 
catdiac  subdivisions  will  be  designated  as  the  standard  Action  Potential 
Distribution  File.  Upon  initialization,  this  vile  is  copied  to  the  current 
Action  Potential  Distribution  File  where  it  can  be  altered  by  the  Action 
Potential  Modification  Program  which  provides  any  change  in  the  current 
profiles  of  action  potential  parameters.  The  regional  anatomy  of  infarct, 
ischemia,  injury,  and  hypertrophy  provide  the  basis  for  establishing  a 
library  of  regional  changes  in  the  action  potential  parameters.  The  Current 
Action  Potential  Distribution  File  is  updated  with  the  changes  and  the 
program  process  continues. 
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4.  Current  Transfer  Impedance  File:  The  current  file  used  during  any 
simulation  procedure  is  a  copy  of  a  file  from  the  library  of  torso  transfer 
impedances.  A  catalog  will  be  maintained  for  choosing  a  specific  heart 
location,  contractile  pattern,  torso  shape  and  conductivities.  This  catalog 
will  uniquely  identify  the  specific  Transfer  Impedance  File  in  the  library. 
Provisions  are  made  in  the  simulation  programs  to  load  the  specific  file  with 
an  entry  from  the  catalog. 
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Discussion  of  the  Flow  Chart 


On  start-up,  the  main  program  requests  the  date  and  time  from  the  oper¬ 
ating  system.  Both  will  be  stamped  on  current  files  and  final  output  data 
forms  for  identification  purposes.  In  addition,  the  IDENT  routine  will 
request  a  unique  identification  number  and  a  short  sentence  that  will  appear 
on  the  output  data  form.  The  purpose  of  this  is  to  describe  and  identify  the 
particular  type  of  input  conditions  or  modifications  that  were  used  in  the 
simulation.  For  example,  the  operator  might  specify  a  2-cm  subendocardial 
infarct  in  segment  5  and  identify  the  output  of  the  simulation  by  the  number 
54321.  The  ECGs,  vectors,  body  surface  maps,  heart  and/or  torso  cross  sec¬ 
tions,  etc,  that  appear  in  hard  copy  and  disk  or  tape  storage  would  have  the 
name  of  the  operator,  date,  time,  ID  number,  and  the  descriptive  sentence  as 
a  label  or  header.  The  main  program  consists  of  a  series  of  six  questions 
for  which  either  yes  or  no  answers  are  given  by  the  operator  at  the  terminal. 

Question  1.  If  the  operator  wishes  to  perform  a  simulation,  he  answers 
yes  to  the  question.  In  the  process  of  initializing  by  responding  to  a 
series  of  questions  in  this  section,  the  current  input  data  files  are  filled 
with  either  standard  or  specific  named  data  files.  Named  data  files  are 
those  files  for  heart  geometry,  Purkinje  network,  action  potential  profiles, 
and  torso  transfer  impedances  that  were  created  in  previous  simulation  ses¬ 
sions  and  saved  on  mass  storage  devices  under  a  specific  name  and  file 
specification.  In  this  way,  a  particular  heart  geometry  with  an  infarct,  for 
example,  could  be  simulated  under  various  other  input  data  changes  and  where 
the  entire  study  might  extend  over  many  sessions  or  days  without  the  need  to 
ever  recreate  the  infarct  geometry  within  the  heart  geometry.  Standard  and 
various  named  input  data  files  would  be  useful  in  debugging  and  testing  the 
simulation  process  following  any  program  modification  or  upgrades. 

Question  2.  If  an  infarct,  hypertrophy,  or  other  heart  geometry  change 
is  desired,  a  yes  response  to  question  2  will  invoke  the  Heart  Geometry 
Alteration  Program.  This  program  requires  the  operator  to  identify  the 
regions  to  be  altered  and  then  input  the  changes  through  the  digitizing  table 
from  a  drawing  or  photograph  of  the  heart  geometry  under  modification,  usu¬ 
ally  the  current  standard  heart  geometry  file  data.  This  alteration  program 
updates  the  current  Heart  Geometry  File  to  reflect  the  modifications  and  then 
returns  to  the  Main  Program. 

Question  3.  A  yes  response  invokes  the  Purkinje  Distribution  Modifica¬ 
tion  Program.  In  addition  to  modifying  the  Purkinje  distribution  over  the 
subendocardial  surface  this  program  is  used  to  specify  the  starting  points 
for  the  Purkinje  activation.  Such  starting  points  correspond  to  the  posi¬ 
tions  within  the  Purkinje  network  where  the  arborizations  of  the  right  and 
left  bundles  of  the  specialized  conduction  system  are  inserted.  Hence,  the 
simulation  of  certain  conduction  defects  is  accomplished  by  first  modifying 
the  start  points  with  this  Purkinje  Modification  Program. 

Changing  the  Purkinje  network  is  accomplished  in  the  same  way  as  the 
heart  alterations,  that  is  by  digitizing  the  additional  regions  or  deleting 
regions  that  are  damaged.  Special  consideration  and  care  must  be  used  when 
modifying  the  Purkinje  network  whenever  the  heart  geometry  is  altered  at  the 
endocardial  surface.  By  definition,  the  Purkinje  network  is  on  the  endocar¬ 
dial  surface  and  must  be  contiguous  with  the  subendocardial  myocardium 


everywhere.  In  addition,  the  network  must  be  continuous  at  a  1  millimeter 
spacing  to  represent  normal  Purkinje  and  yield  a  normal  excitation  velocity 
in  the  wave  propagation  program.  These  considerations  will  be  covered  in 
detail  in  the  user  guide.  After  inputing  the  new  starting  points  or  modifying 
the  Purkinje  the  program  updates  the  current  Purkinje  Distribution  File  and 
returns  to  the  Main  Program. 

Question  4.  A  yes  response  invokes  the  Action  Potential  Distribution 
Modification  Program.  This  program  does  not  alter  the  definition  of  the 
cardiac  segments  and  zones.  The  active  myocardium  is  subdivided  into  20 
segments,  12  in  the  left  ventricle  and  septum  and  8  in  the  right  ventricle. 
Each  segment  is  further  subdivided  into  zones  from  endocardium  to  epicardium. 
Using  drawings  of  the  segments  and  zones  within  the  segments  the  operator 
responds  to  requests  for  specific  regions  or  zones  where  action  potentials 
are  to  be  altered.  The  operator  then  enters  the  new  action  potential  para¬ 
meters  for  the  particular  local  region.  Requests  for  additional  input  con¬ 
tinue  until  the  operator  indicates  the  end.  The  Action  Potential  Distribu¬ 
tion  Program  updates  the  current  Action  Potential  Distribution  Files  and 
returns  to  the  Main  Program. 

Question  5.  Specifying  that  a  new  set  of  transfer  impedances  from  the 
cardiac  generator  to  the  torso  surface  be  introduced  into  the  simulation  is 
accomplished  by  table  look-up  of  the  desired  file  from  the  library  of  trans¬ 
fer  impedances.  The  library  is  created  by  solving  the  appropriate  boundary 
value  problem  for  various  torso  shapes,  lung  and  blood  mass  geometries,  heart 
locations  and  rotations,  regional  wall  motion  abnormalities  and  internal 
conductivities.  A  large  set  of  boundary  value  solutions,  containing  known 
important  differences  in  geometry  and  conductivity  will  be  computed  and 
uniquely  labeled  and  described.  This  catalog  and  disk  library  will  be  updated 
as  additional  sets  of  transfer  impedances  arrays  are  produced  during  the  02 
(and  03)  year(s). 

Question  6.  Under  certain  conditions  such  as  debugging  or  testing  of  the 
program,  it  will  be  desirable  to  enter  the  utilities  program  without  exe¬ 
cuting  the  simulation.  A  no  response  to  question  6  would  invoke  the  util¬ 
ities  program  and  allow  the  operator  to  rename  files,  print  or  plot  data,  or 
other  utility  type  operations  invoked  through  operating  system  requests. 
Otherwise,  a  yes  response  will  send  the  simulation  programs  into  execution  at 
the  appropriate  entry  points  defined  by  the  responses  to  the  previous  ques¬ 
tions.  If  the  operator  enters  the  utilities  program  at  question  6,  it  may  be 
desirable  to  chose  between  exiting  the  program  or  branching  back  to  the 
starting  point  in  the  main  routine.  This  re-entry  question  will  have  to  be 
considered  in  further  detail  at  the  time  of  the  program  generation  because  of 
the  many  ramifications  of  re-entry  as  relates  to  the  operating  system 
characteristics. 


Simulation  Programs 

Positions  3,  4,  and  5  in  the  Main  Program  are  entry  points  in  the  simula¬ 
tion  determined  by  answers  to  Questions  1-5.  The  program  structure  mini¬ 
mizes  the  execution  time  by  not  using  the  entire  set  of  programs  when  no 
changes  have  occurred.  Entry  at  point  3  causes  the  entire  set  of  programs  to 
be  executed  starting  with  the  wave  propagation  program.  Current  data  files 


60 


shown  an  the  left  in  the  flow  diagram,  and  generated  output  files  appear  on 
the  right-hand  side.  As  illustrated  in  the  flow  chart,  the  Wave  Propagation 
Program  inputs  the  current  Heart  Geometry  and  Purkinje  Network  Files  and 
produces  as  output  the  current  Cardiac  Activation  Data  File.  Cardiac  Activa¬ 
tion  Data  File  becomes  the  input  to  the  next  program,  the  Equivalent  Cardiac 
Sources  Program  which  combines  it  with  the  current  Action  Potential  Distribu¬ 
tion  Data  File  and  generates  an  Equivalent  Cardiac  Source  File  representing 
the  electrical  activity  of  the  heart  being  simulated  during  depolarization 
and  repolarization. 

The  last  stage  of  the  simulation  process  is  the  ECG  Program  which  com¬ 
bines  the  Cardiac  Source  and  the  volume  conductor  Transfer  Impedance  Data 
Files  to  produce  the  output  ECG's,  VCG's,  and  body  surface  maps.  The  latter 
may  be  displayed  as  isopotential  or  perspective  plots  for  each  2.5  msec  or  as 
isoarea  maps,  1  plot  each  for  QRS,  ST,  T,  or  QRST. 

The  utility  program  is  inserted  at  the  end  of  the  simulation  program  to 
permit  the  operator  to  exercise  various  options  for  outputting  data  or 
storing  data  for  future  reference.  Some  of  the  routines  within  the  utility 
program  would  be  system  requests  such  as  copying,  renaming,  or  deleting  files 
from  the  mass  storage  device.  Other  options  would  include  printing  ECG  data, 
plotting  ECG  data,  printing  or  plotting  any  of  the  current  input  data  files. 
As  the  forward  model  programs  are  developed,  the  utilities  required  to  test, 
debug,  and  analyze  the  simulation  process  will  be  added  to  the  utility  pro¬ 
gram.  Various  graphics  aids,  given  appropriate  hardware,  could  also  be 
added. 


Color  Heart  Display  Program 

This  program  will  be  developed  under  Task  2A  and  consists  of  a  color 
graphics  display  of  the  heart  geometry  under  simulation  at  the  time  it  is 
called  into  execution.  Two  input  data  files  are  processed  into  the  display: 
the  current  heart  geometry  and  the  current  action  potential  distribution 
files.  The  operator  will  select  from  a  menu  the  slices  of  the  heart  to  be 
displayed.  From  the  above  2  input  files,  the  heart  geometry  will  be  dis¬ 
played  on  a  dark  green  background.  Healthy  myocardium  will  be  presented  in 
red,  infarcted  regions  in  white,  and  ischemic  zones  in  blue. 


Input  Data  Modification  Programs 
and  Enhancements' 


Initialization  Program 

Generally,  the  initialization  program  is  the  first  program  to  be  exe¬ 
cuted.  Its  purpose  is  to  fill  the  current  data  files  with  either  a  standard 
or  normal  configuration.  A  set  of  standard  or  normal  input  data  files  is 
used  to  simulate  normal  ECGs  or  generate  a  standard  test  set  of  output  data. 
These  files  are  required  for  debugging  the  program  following  any  modifica¬ 
tions  to  the  source  code.  In  addition,  after  a  simulation  session  on  the 
computer,  the  current  data  files  can  be  renamed  and  saved  with  the  utility 
program.  Therefore,  the  next  simulation  session  could,  on  initialization. 
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specify  any  of  the  previous  files  to  be  reloaded  into  current  data  files.  At 
the  start  of  a  simulation  session,  the  operator  should  assume  the  current 
data  files  to  be  undefined  and  execute  the  initial ization  routine,  thereby 
creating  new  current  data  files  with  standard  or  normal  data,  or  by 
specifying  some  file  or  files  that  had  been  preserved  from  previous  sessions. 


The  initialization  program  will  be  designed  to  solicit  responses  from  the 
operator  as  to  which  files  to  load  or  default  to  the  normal  simulation  data. 
For  example,  a  normal  set  of  transfer  impedances  from  the  cardiac  generators 
to  the  ECG  leads  at  the  chest  surface  would  always  be  loaded  by  default,  a 
normal  adult-male  heart  geometry  for  the  wave  propagation  program  with  its 


normal  Purkinje  distribution. 


set  of  distributed  action  potential 


parameters.  These  files  are  required  to  exist  at  the  time  of  execution  and 
would  be  set  to  normal  values  unless  decided  otherwise  by  the  operator.  To 
accomplish  this,  the  program  will  request  the  filespec  of  the  desired  current 
data  files  at  the  beginning  of  the  program. 

Geometry  Alteration  Program 

Figure  13  shows  the  form  in  which  the  heart  geometry  is  represented. 
This  particular  cross-section  happens  to  be  at  location  Z=51  mm,  measured 
from  the  apex  towards  the  base  of  the  heart.  A  similar  cross-section  exists 
for  each  1  mm  of  myocardium.  In  the  horizontal  projection  the  heart  geometry 
is  defined  by  a  raster  scan  in  the  Y  direction;  that  is,  for  each  1-mm 
increment  in  the  Y  direction  the  values  of  the  Y  coordinates  for  fixed  X 
specify  the  endocardial  and  epicardial  surfaces.  In  Figure  13,  for  example, 
the  value  Xn  represents  a  fixed  value  of  the  X  coordinate  and  the  scan  in  the 

Y  direction  picks  up  the  values  Y1,Y2, . Y6.  The  myocardium  in  which  the 

wave  of  excitation  will  propagate  is  input  in  the  following  format: 

Z= 5 1 

(Xn,Yl,Y2) 

(Xn,Y3,Y4) 

(Xn ,Y5 ,Y6) . 

With  this  input  data  the  wave  propagation  program  will  construct  "normal 
myocardium"  between  all  coordinates  between  Y1  and  Y2  for  the  given  value  of 
Xn  and  Z=50.  All  other  coordinates  are  set  this  same  way,  leaving 
unspecified  coordinates  as  volume  regions  where  the  wavefront  cannot 
propagate. 


If  infarction  is  to  be  simulated,  then  the  operator  must  specify  the 
infarct  in  the  same  manner  through  the  heart  geometry  alteration  program. 
For  the  purpose  of  the  simulation,  an  infarct  specification  is  defined  as  a 
region  of  the  myocardial  geometry  where  the  activation  front  will  not  be 
allowed  to  propagate.  This  means  that  the  activation  front  will  propagate 
around  the  3-dimensional  infarct  region  and  the  infarct  will  not  contribute 
to  the  equivalent  cardiac  generators  producing  the  body  surface  ECGs, 

As  an  example  of  how  the  first  infarcts  were  simulated,  consider  the 
cross-section  of  normal  heart  geometry  in  Figure  14  taped  to  the  digitizing 
ta^e.  An  infarct  was  defined  by  a  line  circumscribing  the  region.  The  digi¬ 
tizing  cursor  was  reset  at  position  1  to  define  the  origin;  then  cursor 


Figure  13.  Example  of  the  raster  scan  that  defines  the  cross  section  of 
cardiac  geometry. 


position  readings  were  taken  at  locations  2  and  3  to  correct  for  possible 
rotation  of  the  X  and  Y  axes.  Next  the  operator  scanned  the  infarct  region 
recording  the  coordinate  positions  (X1,Y1,Y2),  (X2, Y1 , Y2) , . . . (Xk ,Y1 , Y2) .  A 
request/answer  program  prompted  the  operator  for  the  Z  level  of  the  heart 
geometry  or  to  exit.  Following  the  last  level  digitized  the  program  conca¬ 
tenates  the  infarct  geometry  file  to  the  normal  heart  geometry  file  to 
produce  an  input  file  for  the  wave  propagation  program. 

During  Task  2A  the  simulation  program  will  insert  infarct  data  manually 
by  appending  the  infarct  coordinates  to  the  current  heart  geometry  file  that 
was  loaded  initially  with  the  normal  geometry. 

A  menu-driven  method  of  inserting  infarction  into  the  zones  within  each 
segment  will  be  provided  during  Task  2B.  Those  combinations  of  infarction 
data  the  operator  selects  will  be  appended  to  the  current  heart  geometry  file 
from  a  library  containing  infarct-specification  coordinates  for  every  zone  of 
each  segment.  For  example.  Figure  15  will  be  displayed  when  level  61  is 
selected  by  the  operator  at  the  keyboard.  The  operator  will  be  requested  to 
enter  the  numbers  corresponding  to  the  desired  region  of  infarction. 


Figure  14.  The  geometry  of  a  typical  infarct  is  superimposed  on  the  geometry 
of  an  appropriate  heart  cross-section. 


An  automatic  form  of  interactive  graphics  hardware/software  such  as  a 
light  pen  or  plotter  and  digitizer  will  be  developed  in  Task  3  to  alter  heart 
geometry.  This  will  permit  the  operator  to  display  sections  of  the  heart  on 
a  graphics  terminal,  and  interacti vely  specify  infarction  on  the  displayed 
cross-section  or  modify  the  geometry  to  include  dilation  of  the  ventricles. 
Such  a  program  will  be  designed  to  automatical ly  append  either  infarction 
and/or  hypertrophy  to  whatever  data  is  in  the  current  heart  geometry  file. 
Figure  16  shows  the  same  cross-section  as  before;  however,  the  operator  is 
requested  to  trace  out  any  infarcted  region  and  trace  out  any  hypertrophied 
region  to  be  included  as  illustrated  in  the  drawing.  Obviously,  the  operator 
must  have  available  in  advance  the  information  on  where  to  draw  the  regions. 


urkinje  Modification  Program 


Geometry  for  the  Purkinje  distribution  is  taken  from  the  endocardial 
surface  of  the  normal  human  heart  geometry.  It  includes  all  of  the  endocar¬ 
dial  surface  with  the  exception  of  the  region  near  the  base  of  the  heart. 
There  are  only  several  limited  conditions  under  which  the  Purkinje  network 
can  be  altered  or  modified.  One  condition  for  Purkinje  modification  will 
occur  if  the  endocardial  surface  of  the  heart  geometry  is  altered  by  the 
geometry  modification  program.  If  such  a  heart  geometry  modification  extended 
the  endocardial  surface  either  into  cavity  or  decreased  the  wall  thickness 
from  the  endocardial  surface,  then  that  portion  of  the  Purkinje  distribution 
effected  would  have  to  be  modified  by  this  routine.  The  reason  for  this  is 
that  the  Purkinje  distribution  must  be  contiguous  with  the  endocardial 
surface  geometry. 

A  second  reason  for  modifying  the  Purkinje  network  would  be  to  simulate  a 
sub-endocardial  infarction  for  which  the  Purkinje  network  under  the  infarc¬ 
tion  has  not  survived.  Therefore,  the  activation  of  the  Purkinje  would 
propagate  around  the  infarcted  region. 

Purkinje  input  data  has  the  same  format  as  the  raster  scan  of  the  heart 
geometry.  Modifying  the  Purkinje  input  data  follows  the  same  procedure  as 
for  the  heart  geometry.  A  plot  of  the  normal  heart  geometry  is  mounted  on 
the  digitizing  table.  Regions  of  the  Purkinje  network  to  be  redefined  are 
outlined  and  then  the  cursor  scans  in  a  raster  mode  and  picks  up  the  points. 
The  operator  must  exercise  care  and  be  certain  the  Purkinje  network  is  con¬ 
tinuous  everywhere  at  a  1  mm  resolution;  that  is,  there  can  be  no  distances 
between  points  of  the  Purkinje  network  that  exceed  1  mm. 

Starting  points  are  selected  from  the  Purkinje  distribution  data  and 
inserted  at  the  end  of  the  current  Purkinje  distribution  file.  A  set  of 
normal  starting  points  will  be  used  for  default  values  if  not  specifically 
entered  by  the  operator. 

Task  2A  will  include  a  set  of  starting  points  and  a  normal  Purkinje 
distribution  on  the  endocardial  surface  of  the  normal  heart  geometry,  with 
changes  made  manually  to  the  Purkinje  distribution  file. 

A  menu-driven  program  to  specify  regional  Purkinje  infarction  is  probably 
not  justified  as  an  interim  solution  for  modifying  the  network,  since  little 
could  be  accomplished  with  the  simulation  program  for  such  variations  in 
Purkinje  morphology. 

During  Task  3,  the  same  interactive  graphics  method  of  modifying  the 
heart  geometry  will  be  incorporated  into  the  design  for  modifying  Purkinje 
distribution  and  specifying  starting  points. 


Action  Potential  Distribution 


Variations  in  the  action  potential  distribution  are  simulated  by  modi¬ 
fying  the  current  activation  potential  file  containing  the  distributed  vol¬ 
tage  profiles.  Ischemia  and  other  action  potential  related  anomalies  are 
simulated  by  specifying  the  coefficients  in  the  action  potential  profiles  as 
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they  apply  to  regions  of  the  heart.  A  discussion  of  the  coefficients  and 
equation  for  the  potential  profile  is  contained  in  the  Forward  Model 
Treatise.  All  regions  of  the  non-infarcted  heart  geometry  require  a  set  of 
coefficients  before  the  equivalent  cardiac  generators  can  be  calculated. 
During  Task  2A  of  the  program  development,  a  set  of  normal  coefficients  for 
every  zone  of  each  segment  of  the  heart  subdivision  will  be  established  and 
incorporated  into  the  model  as  the  default  values.  Modification  of  the 
regional  values  of  the  coefficients  will  be  menu-driven  for  the  various  zones 
of  each  segment  as  depicted  in  Figure  15,  a  cross  section  of  the  heart 
showing  the  configuration  of  a  mid-ventricular  slice.  The  program  will 
prompt  the  operator  for  the  segment  number  followed  by  a  prompt  for  the  zone 
within  that  segment.  The  operator  will  insert  the  value  of  the  coefficient 
requested  by  the  program.  When  no  further  coefficients  are  to  be  specified 
the  operator  exits  the  routine,  which  inserts  the  modification  into  the 
current  action  potential  distribution  file  for  input  to  the  equivalent 
cardiac  source  program. 

An  automatic  method  of  specifying  the  regions  of  ischemia  to  be  simulated 
is  scheduled  for  development  during  Task  2B.  The  proposed  method  will 
display/plot  a  cross-section  of  the  heart  geometry,  contained  in  the  current 
file,  where  the  ischemic  regions  are  to  be  specified.  The  operator  will  label 
the  region  uniquely,  then  with  a  light  pen  or  cursor  trace  out  the  region. 
At  the  conclusion  of  tracing,  the  program  will  prompt  the  operator  for  a  set 
of  coefficients  to  be  assigned  to  the  action  potential  profile  for  this 
unique  region.  The  program  will  loop  in  this  mode  until  the  operator  has 
finished  all  modifications.  Again,  this  data  will  be  inserted  into  the 
current  action  potential  distribution  program.  Figure  17  depicts  the  input 
design  for  interactive  specification  of  ischemic  geometry. 


Transfer  Impedance  Modification  Program 


The  proposed  transfer  impedance  modification  rou 
look-up  for  specific  cases.  A  set  of  transfer  impeda 
as  the  solution  for  a  boundary  value  problem  for  any 
try,  cardiac  source  locations,  and  conductivities, 
impedances  will  be  computed  throughout  Tasks  2A,  2B, 
impedance  selection  routine  will  prompt  the  operator 
tion  of  the  desired  set  of  transfer  impedances.  Fil 
in  a  catalog  of  the  updated  library  and  include  a  un 
transfer  impedances. 
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Simulation  of  wall  motion  during  systole  and  cardiac  dilatation  are 
examples  that  require  specification  of  the  cardiac  source  locations  as  a 
function  of  time.  A  discussion  of  how  this  information  must  be  specified  is 
in  the  ECG  Forward  Model  Treatise.  Briefly,  a  set  of  cardiac  generator 
locations  must  be  specified  at  several  instants  of  time  during  the  systolic 
interval,  the  solution  to  the  boundary  value  problem  is  computed  for  each 
instant,  and  then  the  remaining  transfer  impedances  are  computed  by 
interpolating  over  this  finite  set. 


Figure  17.  Same  cross-section  and  infarct  is  shown  as  Figure  16.  A  border 
zone  of  ischemia  is  shown. 


Enhancements  to  Forward  Model 

Several  enhancements  to  the  wave  propagation  simulation  program  are  under 
consideration  as  part  of  Task  3: 

1.  Rather  than  execute  the  program  in  two  separate  phases,  first 
the  rapid  Purkinje  activation  followed  by  slower  cell-to-cell  conduction, 
the  program  can  be  modified  to  run  both  Purkinje  and  wall  activation 
concurrently  at  their  respective  velocities.  Such  an  enhancement  would 
permit  the  simulation  of  conduction  system  defects  that  result  in 
Purkinje  re-entry.  This  enhancement  will  be  added  to  the  original 
activation  model . 

2.  Another  enhancement  to  the  propagation  model  would  be  to 
include  non-uniform  conduction  velocity  according  to  certain  regional 
specifications,  because  ischemic  regions  conduct  slower  than  normal  myo¬ 
cardium  and  activation  in  normal  tissue  is  asymmetric.  Two  concerns 
mediate  against  such  an  enhancement  (over  and  above  the  obvious  increase 
in  computer  resources):  a)  the  information  from  laboratory  experiments  is 
not  available  to  specify  the  magnitude  and  direction  of  the  conduction 
velocity  in  the  intact  human  heart;  and  b)  it  hasn't  been  ascertained 
whether  the  inclusion  of  non-uniform  velocity  is  a  first-order  effect 


that  projects  onto  the  body  surface  EC6  leads.  Experiments  will  be 
conducted  to  determine  if  the  simulation  of  regional  variations  in  con¬ 
duction  velocity  produce  a  first-order  effect  in  the  surface  ECGs.  If 
the  results  demonstrate  that  velocity  variations  are  first-order,  then 
the  propagation  model  will  be  enhanced  to  include  both  variations  in 
conduction  velocity  and  Purkinje  re-entry  into  a  single  module. 

3.  During  Task  3,  an  analysis  should  be  performed  to  determine 
whether  there  exists  a  set  of  torso  geometries,  lung  fields,  blood 
masses,  conductivity  ratios  and  generator  locations  for  which  transfer 
impedances  could  be  calculated  and  then  by  interpolation  and/or  extra¬ 
polation,  all  other  transfer  impedances  could  be  obtained  from  a  specifi¬ 
cation  of  parameters,  such  as  A-P  and  lateral  dimensions,  height,  or 
information  from  chest  x-rays,  echos,  and  ventriculograms. 


Numerical  Experiments  with  the  Model 
and  Sequence  of  its  Development 


Systolic  Wall  Motion  and  Dipole  Location 

The  most  computer-intensive  portion  of  the  modelling  effort  is  the  gener¬ 
ation  of  a  library  of  transfer  impedances  for  the  various  combinations  of 
torso  resistivities  and  heart  locations.  This  is  due  to  the  fact  that  the 
Gelernter  Swihart  (G-S)  simulation  of  the  homogeneous  torso  volume  conductor 
requires  the  solution  of  some  3000  simultaneous  equations.  It  therefore 
becomes  of  primary  importance  to  investigate  the  relative  importance  of 
changes  in  heart  location  that  might  be  related  to  systolic  wall  motion  and 
the  potential  significance  of  exercise-related  changes  in  the  volume  conduc¬ 
tor  early  on  in  the  systems  analysis.  Given  that,  at  the  source,  one  would 
expect  a  10%-20%  change  in  ST  segments  opposite  to  the  major  QRS  deflect, on 
as  a  reasonable  indicator  of  subendocardial  ischemia,  a  series  of  experiments 
can  be  performed  in  the  first  months  of  this  effort  to  begin  to  assess  if  the 
kind  of  regional  wall  motion  changes  expected  with  systole,  including  the 
rotational  effects,  are  first-order  or  lesser-order  effects. 

Specifically,  the  potential  distribution  on  a  typical  adult  male  torso  (a 
body  surface  ECG  map)  from  the  propagation  model  for  QRS  will  be  computed. 
Antero-septal ,  superoapical ,  posteroapical ,  inferobasal  and  posterobasal  unit 
dipoles  from  this  model  will  be  moved  10  mm  in,  out,  basally,  apical ly, 
clockwise,  and  counter-clockwise  from  the  expected  median  position;  G-S 
solutions  to  the  inhomogeneous  volume  conductor  will  be  calculated  for  each 
perturbation  in  location;  the  absolute  change  in  voltage  that  would  be  pro¬ 
duced  on  a  surface  map  by  a  10%-20%  change  in  each  unit  dipole  will  be 
assessed.  In  the  analysis  of  the  data  it  will  be  asked:  Which  dipoles  when 
changed  by  either  10%  or  20%  provoked  changes  in  surface  measured  voltages 
outside  the  expected  noise  in  clinically  recorded  rest  and  exercise  ECG 
surface  maps?  A  relevant  cardiological  and  clinical  question  is:  What  is  the 
extent  of  regional  subendocardial  ischemia  that  would  be  indicative  of  sig¬ 
nificant  coronary  narrowing  in  each  of  the  major  coronary  distributions?  For 
this  purpose,  the  main  diagonal  branch  of  the  left  anterior  descending  coro¬ 
nary  artery  will  be  considered  as  a  separate  arterial  distribution  since  it 
perfuses  a  volume  of  myocardium  roughly  equivalent  to  that  perfused  by  the 
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right  and  circumflex  coronary  arteries.  Analysis  of  this  data  in  the  context 
of  clinically  relevant  areas  of  ischemia  will  give  a  first  order  approxima¬ 
tion  of  the  potential  resolution  of  this  method.  Earlier  outputs  from  the 
model  (see  Figure  11  in  the  Background  section)  suggest  that  significant 
subendocardial  ischemia  involving  the  major  portion  of  the  inner  2-3  mm  of 
each  of  the  12  LV  segments  will  produce  changes  above  the  noise  level  that 
can  be  achieved  in  the  typical  exercising  patient  with  normal  torso  geometry 
and  resistivities. 

A  related  question  that  might  result  in  very  major  savings  in  computation 
time  is:  What  are  the  potential  physiological,  and  clinical  cardiological 
results  of  considering  these  12  regions  as  having  "fuzzy  borders"?  What  are 
the  conceptual  electrophysiological  consequences  of  assuming  the  borders  have 
that  same  10-mm  uncertainty  by  keeping  the  unit  dipole  location  constant  and 
not  trying  to  simulate  rotational  and  motion  changes  with  systole,  and  with 
the  dilatation  of  ischemia?  It  may  be  quite  enough  to  say  that  there  is 
significant  ischemia  somewhere  in  each  of  the  4  main  coronary  artery  distri¬ 
butions.  The  main  question  now  would  be:  When  this  simplifying  assumption  is 
made  and  a  single  set  of  transfer  matrices  is  used,  does  the  failure  to  be 
realistic  at  the  source  level  result  in  interpretable  body  surface  maps,  and 
can  the  validation  criteria  outlined  on  pages  52  and  53  be  met? 


Subendocardial  Ischemic  Area  Size: 
Sensitivity  and  Resolution  Studies 


When  the  action  potential  distribution  profiles  have  been  incorporated  in 
the  whole  ventricular  geometry,  it  will  be  possible  to  approach  this  same 
problem  with  a  more  rational  and  physiologically  quantitative  approach.  For 
example,  one  can  now  ask:  How  would  subendocardial  ischemia  involving  the 
inner  2  mm  of  each  region  but  only  extending  over  an  elliptical  area  mea¬ 
suring  10  x  15  mm  be  seen  on  the  body  surface?  This  could  be  assessed  with  a 
single  typical  set  of  transfer  impedances  for  the  mid  systolic  position.  The 
size  of  the  internal  surface  area  of  ischemia  and  its  transmural  extent  could 
be  systematically  varied  to  assess  the  limits  of  resolution  of  each  region  on 
the  body  surface,  and  which  local  ECG  leads  best  discriminate  that  change. 
What  size  or  extent  of  typical  subendocardial  ischemia  in  each  of  the  4  main 
coronary  artery  distributions  can  be  identified  in  the  12  lead  ECG  when  the 
simulation  of  the  non-ischemic  normal  is  available  for  comparison,  and  what 
criteria  in  these  leads  would  identify  this  regional  ischemia?  A  typical 
experiment  simulating  subendocardial  ischemia  in  the  distribution  of  the 
distal  circumflex  and  the  marginal  circumflex  coronary  arteries  is  shown  in 
Figure  18. 


Intramural  Conduction  Velocities  with  Ischemia 


While  the  data  is  conflicting  about  whether  during  acute  subendocardial 
ischemia  the  conduction  velocity  through  the  wall  is  slowed  or  accelerated, 
both  effects  can  be  explored  systematically  with  the  model.  To  be  efficient 
about  these  numerical  experiments,  these  will  best  be  accomplished  after  the 
"Input  Data  Modification  Programs  and  Enhancements"  described  on  pages  61-69 
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Figure  18.  A  typical  numerical  experiment  evaluating  the  effect  of  local  sub¬ 
endocardial  ischemia  in  the  region  of  the  distal  circumflex  coro¬ 
nary  artery  (segment  12)  and  the  marginal  circumflex  (segment  11) 
is  shown  in  this  diagram.  The  effects  of  either  alone  on  the  body 
surface  potential  distribution  as  compared  to  normal  (ECG  surface 
map  changes)  can  be  assessed.  The  effect  of  proximal  circumflex 
coronary  narrowing  can  be  simulated  by  combining  the  effects  of 
both  of  these  two  segments.  Regional  wall  motion  abnormalities  in 
these  same  regions  can  be  simulated  by  using  appropriate  transfer 
matrices  from  the  Transfer  Matrix  Library  that  emulate  dilatation 
and  regional  thinning  in  the  local  segment. 

have  been  accomplished.  To  aid  in  the  planning  and  priority  setting  of  the 
development  of  the  model,  preliminary  assessment  of  the  effect  of  varying 
local  intramural  conduction  velocities  will  be  done  late  in  the  1st  year  of 
this  development  effort. 


Combined  Effects  of  Wall  Motion  Abnormalities  and  Ischemia 

When  the  studies  described  above  have  been  completed,  including  simu¬ 
lating  large  multi-segment  areas  of  subendocardial  ischemia,  the  next  phase 
will  be  to  combine  location  changes  with  regional  ischemic  changes.  The  plan 
is  to  change  each  variable  systematical ly.  We  will  first  assess  the  role  of 
regional  thinning  and  dilatation  with  ischemia,  which  has  been  shown  to  be 


one  of  the  earliest  changes.  We  will  then  compare  what  this  change  does  to 
the  ECG  surface  map  with  that  produced  by  action  potential  changes  in  sub¬ 
endocardial  regions,  as  separate  effects  and  as  combined  effects  (for  small, 
medium,  large,  and  very  large  areas  of  regional  ischemia).  Most  of  these 
systematic  studies  done  by  varying  generator  locations  will  depend  on  having 
appropriate  sets  of  transfer  impedances  in  the  library.  These  will  take  time 
to  accumulate.  To  make  the  systematic  modifications  in  the  heart  geometry 
files,  and  the  action  potential  distribution  files  in  an  efficient  manner 
will  require  that  the  input  data  modification  programs  be  available.  These 
studies  will  therefore  be  done  in  the  main  after  the  development  of  the 
computer  routines  to  modify  the  input  data  to  each  of  the  stages  of  the 
program.  This  is  Subtask  28,  page  51  (discussed  in  more  detail  on  pages  61- 
69)  and  is  projected  to  be  done  in  the  2nd  year  of  this  development  effort. 


Combined  Infarction,  Injury,  and  Ischemia 

With  the  comprehensive  set  of  programs  for  easily  modifying  the  input 
data  files  we  will  also  be  in  a  position  to  systematically  vary  the  heart 
geometry  to  simulate  all  sizes  of  infarction  in  any  distribution  from  sub¬ 
endocardium  to  intramural,  to  subepicardial  and  transmural  in  any  region  of 
the  heart  including  the  right  ventricle,  with  or  without  hypertrophy  and/or 
dilatation.  All  variations  of  subendocardial,  intramural,  subepicardial,  and 
transmural  injury  and  ischemia  in  the  various  coronary  artery  distributions 
can  be  simulated  by  using  the  Action  Potential  Distribution  Modification 
Program.  Changes  that  occur  in  any  or  all  of  the  above  in  the  presence  of 
interventricular  conduction  defects  can  be  simulated  with  the  Purkinje  Modi¬ 
fication  Program.  With  time,  when  a  comprehensive  library  of  Torso  Impedance 
Files  is  generated,  the  simulation  can  be  made  to  deal  with  all  combinations 
of  body  builds,  internal  changes  in  resistivity  as  might  occur  with  +G 
loading,  exercise,  and  in  chronic  emphysema. 


ENGINEERING  CONSIDERATIONS:  PORTABLE  MULTILEAD  DIGITAL  ECG  CART 


Present  Major  Commercial  Vendors 


Both  vendor  representatives,  Steve  Koerper  (and  now  Francis  Charbonnier) 
for  Hewlett  Packard,  and  Timothy  Mickelsen  for  Marquette  Electronics  indicate 
a  continued  interest  by  their  companies  in  maintaining  close  contact  with  the 
development  of  any  new  ECG  Cart  System  that  might  result  from  the  current 
effort.  Both  have  indicated  that  they  would  give  serious  consideration  to  a 
request  for  bid  on  such  a  system  when  and  if  such  were  to  come  from  the  Air 
Force  and  the  Department  of  Defense.  The  current  Case  II  Exercise  System  by 
Marquette  is  upgradable  now  to  18  leads  sampled  at  4  kHz  by  a  simple  card 
insert.  Up  to  32  leads  in  a  portable  cart  system  can  be  accomplished  with 
only  modest  engineering  effort.  The  current  Hewlett  Packard  system  is  a 
serial  3  channel  system,  and  can  take  any  number  of  leads  in  sets  of  3  by 
developing  a  switching  box  for  the  lead  array  chosen.  The  ECG  R  &  D  section 
of  the  company  will  soon  be  headed  up  by  Francis  Charbonnier,  who  attended 
the  final  peer  review  session  at  USAFSAM.  He  indicates  that  the  company  will 
consider  carefully  the  potential  marketabi 1 i ty  of  a  larger  simultaneous 
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multi-lead  system  as  this  effort  progresses.  Both  vendors  have  extensive 
service  networks.  The  specific  needs  of  the  USAF  for  ECG  cart  service  in  the 
field  will  need  to  be  looked  at  as  this  investigation  progresses. 


EDL  Laboratories 


Roland  Wyatt,  who  worked  for  years  on  the  engineering  aspects  of  the 
University  of  Utah  Cardiology  Sections  multilead  mapping  system,  has  recently 
moved  to  the  Endotech  Development  Laboratory  at  the  Research  Park  of  the 
University  of  Utah  and  is  marketing  a  32-channel  ECG  recording  system.  Ver¬ 
sions  of  the  system  can  be  fabricated  in  modules  of  32-256  channels.  A  192- 
channel  sample  and  hold  version  with  a  1  kHz  sampling  rate  is  currently  being 
fabricated.  Wyatt  indicates  that  his  company  would  be  quite  responsive  and 
anxious  to  bid  on  a  portable  ECG  cart  system  when  and  if  the  investigative 
effort  of  this  project  identified  the  software  and  hardware  engineering 
requirements  of  such  a  system.  The  EDL  company  has  a  working  relationship 
for  the  servicing  and  repair  of  its  solid  state  devices  with  Stortz  and  Wolf 
who  maintain  a  worldwide  24  hr  solid  state  and  video  service.  Wyatt  states 
that  they  have  made  considerable  use  of  this  service  and  consider  it  quite 
responsive  and  adequate  to  the  customer  needs.  He  considers  the  servicing  of 
a  multilead  ECG  portable  cart  exercise  system  to  be  of  similar  complexity  to 
their  existing  devices,  and  believes  the  servicing  of  the  proposed  system  to 
pose  no  special  problems. 


University-Affiliated  Multilead  ECG 
Systems  Development 

Robert  Guardo  of  the  Biomedical  Engineering  Department  of  the  University 
of  Montreal  has  developed  a  highly  buffered  system  in  multiples  of  64  chan¬ 
nels  that  can  be  selected  for  frequencies  from  DC  at  the  low  end  to  4  kHz  at 
the  high  end  of  the  spectrum.  His  section  has  produced  a  limited  number  of 
these  units  for  other  investigators,  but  is  not  planning  to  go  into  a  high 
level  production  mode.  They  would  be  interested  in  a  development  relation¬ 
ship  with  one  of  the  commercial  vendors,  however,  in  the  event  that  a  speci¬ 
fic  need  for  a  large  number  of  portable  carts  evolves  from  the  present 
USAFSAM/Rancho-USC  modelling  project.  The  University  of  Utah/Lux  mapping 
system  currently  undergoing  evaluation  in  the  resting  subject  prior  to  coro¬ 
nary  angiography  at  the  USAFSAM  is  also  adaptable  to  exercise  mapping. 
Neither  this  system  nor  the  University  of  Montreal  one  currently  contains 
software  logic  for  interpretation  of  the  records,  and  this  might  be  the  major 
mutual  benefit  of  either  of  these  groups  developing  a  collaborative 
arrangement  with  industry. 
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INTRODUCTION 


The  forward  model  in  electrocardiography,  presented  in  this  report,  is  a 
mathematical-physical  description  of  the  genesis  of  the  electrocardiogram 
transformed  into  algorithms  and  coded  into  a  computer  program  capable  of 
simulating  a  voltage  distribution  on  the  chest  produced  by  electrical  activ¬ 
ity  in  the  heart.  Basic  laws  in  physics  yield  algorithms  for  determining  how 
the  chest  voltages  are  related  to  the  underlying  electrical  events  in  the 
heart.  The  computer  program  is  the  mechanism  for  expressing  the  algorithms; 
executing  the  program  applies  the  algorithms  to  specific  input  data.  Elec¬ 
trical  events  in  the  heart  and  the  influence  of  the  inhomogeneous  volume 
conductor  are  represented  as  numbers  and  symbols  in  the  computer  program. 
When  the  program  is  executed  on  this  input  data,  the  numbers  and  symbols  are 
processed  in  a  manner  specified  by  the  basic  laws  of  physics,  as  expressed  by 
the  algorithms.  Such  a  computer  program  allows  the  genesis  of  the 
electrocardiogram  to  be  deduced  from  basic  principles. 

Simulating  the  electrocardiogram  with  a  computer  program  is  similar  to 
performing  a  controlled  experiment;  input  conditions  are  specified  to  the 
program  that  correspond  to  the  desired  experimental  subject  and  the  program 
computes  the  resulting  ECG.  Unlike  an  actual  experiment,  where  the  control 
conditions  are  difficult  or  impossible  to  establish,  the  numbers  and  symbols 
in  the  program  are  not  restricted  by  experimental  procedures.  Instead  the 
computer  program  permits  any  input  conditions  that  are  consistent  with  the 
algorithms  and  program  design.  A  forward  model  for  simulating  the  ECG  of 
humans  extends  the  cardiolog’st's  understanding  of  how  the  voltage  distribu¬ 
tion  on  the  chest  is  related  to  normal  and  abnormal  electrical  activity  in 
the  heart  as  well  as  the  influence  of  complicated  torso  volume  conductors. 
Simulation  by  a  computer  may  be  the  only  method  available  to  predict  compli¬ 
cated  ECGs,  develop  hypotheses  about  their  origin,  and  establish  diagnostic 
criteria,  especially  in  those  cases  where  the  abnormal  processes  cannot  be 
identified  apriori  because  of  insufficient  information  on  how  they  are  man- 
fested  in  the  ECG.  Probably  the  most  useful  aspect  of  the  forward  simulation 
program  is  in  its  ability  to  generate  hypotheses,  in  the  form  of  ECGs,  about 
such  abnormal  processes  that  can  subsequently  be  tested  on  a  human 
population. 


DESCRIPTION 

Ohm's  law  for  current  conduction  in  an  extended  medium  is  an  uncompli¬ 
cated  way  to  perceive  the  heart  current  source  imbedded  in  the  inhomogeneous 
thorax  and  the  forward  model.  Stated  as  a  matrix  equation,  a  concise  way  to 
represent  a  large  system  of  equations.  Ohm's  law  becomes: 

V  =  A  x  D  ( A-l ) 

where  V  is  the  matrix  of  body  surface  voltage  distributions,  A  is  a  matrix  of 
impedances  for  the  inhomogeneous  volume  conductor  (commonly  known  as  transfer 
coefficients),  and  D  is  a  matrix  of  current  sources  within  the  heart.  ECGs 
are  contained  as  elements  in  V,  whereas  the  elements  of  A  account  for  the 
effects  of  right  and  left  lung  fields,  blood  masses,  and  the  torso-air 
interface.  Elements  of  D  correspond  to  the  currents  developed  during  depo¬ 
larization  and  repel ari zation  of  myocardial  cells.  Algorithms  in  the 
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simulation  programs  contain  numbers  and  symbols  which  yield  the  elements  of 
these  matrices  when  executed.  The  forward  model  programs  separate  naturally 
into  three  groups:  the  first  group  of  programs  compute  the  elements  of  D, 
secondly,  a  group  of  programs  to  compute  the  elements  of  A,  and  third  a  group 
to  combine  A  and  D  to  yield  the  body  surface  voltages  V.  Computing  the 
elements  of  D  is  unrelated  to  the  solution  of  the  boundary  value  problem  to 
obtain  A.  The  only  restrictions  on  D  is  that  of  a  current  source  obeying  the 
conservation  of  current  law  with  a  specified  location.  Part  1  below  dis¬ 
cusses  the  theory  and  observations  about  D  leading  to  the  algorithms  embodied 
in  the  forward  model.  Part  2  discusses  the  theory  and  approximate  solution 
to  the  inhomogeneous  boundary  value  problem  in  electrocardiography.  Part  3 
brings  together  the  solutions  for  A  and  D  in  a  manner  useful  to  the  cardio¬ 
logist,  i.e.,  in  the  form  of  conventional  12  lead  ECG,  VCGs,  body  surface 
maps,  and  graphic  displays  of  ventricular  activation. 

Programs  for  computing  the  elements  of  D  are  written  mostly  in  Fortran 
with  a  few  assembly  language  subroutines  linked  through  Fortran  calls.  These 
programs  require  a  considerable  amount  of  RAM  memory  to  hold  the  heart  geome¬ 
try,  thereby  avoiding  an  I/O  bound  leading  to  excessive  CPU  overhead.  Exam¬ 
ple:  a  normal -adul t,  male  human  heart  requires  125K  bytes  of  storage  added 
to  a  program  of  130K  bytes.  All  of  the  programs  for  the  boundary  value 
solution  are  written  in  Fortran.  These  programs  fit  easily  into  a  64K  byte 
partition.  Since  the  transfer  impedances  are  an  iterative  solution  to  a  set 
of  2500  simultaneous  equations,  the  programs  are  computationally  intensive. 
Depending  on  the  volume  conductor  specifications  and  stopping  criteria,  a  set 
of  solutions  corresponding  to  various  current  source  locations  will  require 
many  hours  of  CPU  time. 


BASIC  PHYSIOLOGICAL  DEVELOPMENT 


Activation  Sequence 

One  of  the  major  factors  influencing  the  potential  distribution  on  the 
torso  surface  is  the  myocardial  activation  sequence.  Genesis  of  QRS  can  be 
inferred  directly  from  detailed  activation  maps  reconstructed  from  measure¬ 
ments  taken  with  multi-point  intramyocardial  bipolar  electrodes.  Discussions 
of  the  many  temporal  and  spatial  aspects  of  the  activation  sequence  in  the 
dog  and  human  can  be  found  in  the  reports  of  Scher  and  Young  (A-l)  and  Durrer 
(A-2).  It  is  clear  that  much  of  what  is  seen  in  surface  ECGs  and  VCGs  is 
directly  proportional  to  the  area,  direction,  and  timing  of  the  activation 
front. 

Experimental  evidence  for  the  existence  of  the  spread  of  activation  in 
the  human  heart  is  documented  by  studies  of  epicardial  and  intramural  excita¬ 
tion.  Detailed  information  about  human  epicardial  excitation  was  obtained 
from  two  sources:  revived  normal  hearts  and  isochronous  time  relations 
during  surgery.  Small  multi-point  intramural  electrodes  were  used  to  explore 
the  excitation  process  from  the  endocardium  to  the  epicardium.  These  studies 
indicate  minimal  intramural  Purkinje  fiber  penetration  into  the  myocardium. 


Numerous  intramural  and  epicardial  electrogram  recording  methods  have 
been  developed  and  the  electrical  activity  analyzed  in  a  manner  suitable  for 


mathematical  modeling,  especially  the  methods  of  Vander  Ark  and  Reynolds  (A- 
3)  that  observe  the  excitation  front  as  uniform  and  closed  at  a  macroscopic 
level  of  observation.  Micro-electrode  recordings,  both  in  vitro  and  in  vivo 
at  the  epicardium,  reveal  information  on  cellular  potentials  and  nonuniform 
phase  velocity.  However,  no  unifying  derivation  exists  for  expressing  the 
macroscopic  current-vol tage  relationships  in  terms  of  micro-electrode  pheno¬ 
mena.  Looked  at  from  the  perspective  of  a  forward  ECG  simulation,  only  the 
macroscopic  characterization  of  the  electrical  events  within  the  heart  is 
useful.  It  would  not  be  an  efficient  use  of  computer  resources  to  include 
detailed  information  on  the  electrical  activity  that  is  not  manifested  in  the 
torso-surface  potential  distribution.  The  objective  of  the  forward  ECG  model 
is  to  produce  an  accurate  approximation  to  the  biophysical  process  rather 
than  realism  through  inclusion  of  minutiae.  Therefore,  in  vivo  macroscopic 
observations  of  the  electrical  activity  in  the  heart  must  be  recorded  at  a 
level  commensurate  with  the  resolution  of  the  simulation.  It  is  our  choice 
to  take  a  spatial  resolution  of  1  mm  for  the  intramural-epicardial  recordings 
and  the  myocardial  geometry  in  the  simulation  program.  Given  that  the  most 
useful  recording  methods  will  be  those  employing  bipolar  electrodes  having  a 
separation  of  1  mm,  a  series  of  experiments  were  performed  on  baboons  to 
characterize  the  electrical  events  in  the  heart  and  used  as  design  input  for 
the  forward  model . 

Activation  sequence  mapping  with  bipolar  electrodes  is  a  procedure  that 
records  the  passing  of  electrical  excitation  (depolarization)  between  the 
bipolar  electrodes.  Measuring  time  at  which  the  electrical  activity  passes 
over  the  electrodes  and  combining  it  with  the  electrode's  position  in  the 
heart  defines  the  propagation  of  the  depolarization  within  the  myocardium. 
It  is  the  observed  propagation  within  the  myocardium  that  serves  as  a  basis 
for  the  design  of  the  wave  propagation  model  to  simulate  QRS  excitation. 
Because  the  design  and  validation  of  a  propagation  model  require  accurate  and 
extensive  mapping  of  intramural  and  epicardial  activation  sequences,  it 
became  necessary  to  design  and  build  reliable  and  efficient  electrode  arrays 
capable  of  rapid  and  accurate  mapping  of  the  entire  heart.  Two  types  of 
electrode  arrays  were  fabricated:  a  13-point  intramural  needle  and  a  24- 
point  epicardial  array.  The  geometrical  arrangement  of  the  bipolar  elec¬ 
trodes  was  chosen  to  accommodate  the  equations  used  to  analyze  the  recor¬ 
dings.  Figure  A-l  contains  photographs  of  the  two  types  of  arrays:  top 
panel  is  the  intramural  array  and  the  lower  panel  is  the  epicardial  clock 
array. 


Experimental  Observations 


Intramural  multipoint  electrodes  were  designed  to  insure  a  strong  and 
reliable  intramural  bipolar  needle  electrode  configuration;  a  20-gauge  stain¬ 
less  steel  tube  was  used  as  a  cylindrical  support  for  the  wire  electrodes. 
This  20-gauge  tube  had  an  .008-inch  slot  cut  the  length  of  the  tube  which  was 
approximately  3  cm  long.  Opposite  the  slot  in  the  tubing,  13  holes  were 
drilled  through  the  wall  of  the  tubing.  Each  hole  was  .008  inch  in  diameter 
and  spaced  at  1-mm  increments.  Stainless  steel  wire  (304  type)  having  a 
diameter  of  .005  inches  and  coated  with  Teflon  to  a  diameter  of  .007  inches, 
was  threaded  into  the  holes  drilled  through  the  tubing.  The  needle  assembly 
is  sharpened  and  coated  with  epoxy  compound.  The  top  panel  of  Figure  A-l  is 
an  enlarged  view  of  a  complete  needle  electrode. 
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Myocardial  activation  sequences  were  reconstructed  from  the  time  of 
occurrence  of  the  peak  voltage  occurring  between  the  adjacent  (1  mm)  elec¬ 
trodes  and  the  location  of  the  needle  in  the  myocardium.  Activation  sequence 
data  consists  of  50  sites  on  the  heart  surface  where  the  13-point  intramural 
needles  are  inserted.  A  long  intramural  needle  (20  cm)  was  used  to  map  the 
septal  regions  of  the  heart  and  to  obtain  a  limited  amount  of  intramural 
activation  data  prior  to  opening  the  animal's  chest  by  percutaneously 
advancing  the  long  needle  electrode  into  the  myocardium. 

The  clock  electrode  consists  of  12  bipolar  electrodes  located  around  the 
circumference  of  a  1-cm-diameter  circle.  Figure  A-2  indicates  the  12  bipolar 
configuration  and  the  labeling  convention  used  to  designate  the  orientation 
of  the  array  and  position  of  each  bipolar  pair.  Each  recording  pair  consists 
of  .2-mm-diameter,  Teflon-coated  type  304  stainless  steel  wire.  The  exposed 
ends  of  the  wires  are  separated  a  distance  of  1  mm  symmetrically  across  the 
face  of  a  1-cm-diameter  circle.  A  pair  of  wires  are  positioned  every  30 
degrees  around  the  circle.  Epoxy  is  used  to  encase  the  wires  in  place  and 
support  a  handle  made  of  a  soft  metal  tube  such  as  copper  or  aluminum.  The 
24  wires  are  placed  inside  the  tube  for  protection.  A  typical  clock  elec¬ 
tive  designed  for  free-hand  epicardial  mapping  is  shown  in  the  bottom  panel 
of  Figure  A-l.  A  12  o'clock  position  is  marked  on  the  epoxy  support  to 
indicate  the  orientation  of  the  array.  Electrode  pairs  are  numbered  as  the 
hours  on  a  clock  face,  looking  down  onto  the  array  in  its  normal  resting 
position  on  the  epicardial  surface,  as  depicted  in  Figure  A-2. 

A  convention  was  selected  which  assigns  the  positive  or  upright  voltage 
deflection  to  activation  propagating  outward  along  a  radius  from  the  center 
of  the  circle  containing  the  bipolar  electrodes.  Electronic  amplifiers 
measure  the  voltage  at  the  outside  electrode  wire  with  respect  to  the  inside 
wire. 


Figure  A-3  is  a  recording  of  the  electrograms  generated  across  the  face 
of  the  clock  electrode  on  the  epicardial  surface  of  a  normal  baboon  heart. 
This  recording  is  from  a  site  on  the  left  ventricle  near  the  apex.  For  the 
purpose  of  describing  the  basis  underlying  the  propagation  model  of  cardiac 
activation,  the  baboon  is  presented  as  a  representative  application  rather 
than  a  study  of  the  activation  sequence  in  the  baboon. 

In  addition  to  the  12  bipolar  recordings,  two  reference  signals  were 
recorded  from  sites  selected  on  the  free  wall  of  the  right  ventricle  near  the 
pulmonary  outflow  tract.  Some  form  of  reference  signal  is  required  to  time- 
align  the  electrograms  recorded  from  various  sites  on  the  epicardial  surface. 
Various  reference  recording  sites  were  tested,  such  as  torso  surface  ECG, 
intramural  bipolar  and  intramural  unipolar  leads.  Our  experience  shows  that 
unipolar  and  bipolar  leads  placed  at  the  outflow  tract  of  the  right  ventricle 
provided  the  most  stable  and  trouble-free  sites  for  recording  reference 
electrograms.  This  region  of  the  heart  produced  repeatable  complexes  while 
the  physical  position  of  the  heart  was  altered  slightly  during  the  mapping 
process.  In  addition,  this  location  minimized  the  accidental  dislocation  of 
the  reference  electrodes.  Separation  on  the  bipolar  reference  electrode  was 
held  at  1  mm. 

At  each  bipolar  lead  of  the  clock  electrode  the  arrival  time  of  the 
activation  front  at  the  center  of  the  bipolar  pair  was  chosen  at  the  peak 


value  of  the  electrogram.  For  each  of  the  12  pairs,  there  is  a  time  tn, 
where  n  =  1,2, ,,,,12  corresponding  to  the  arrival  of  the  activation  front  at 
the  electrode  pair.  This  time  t  n  is  calculated  from  a  characteristic  point 
on  either  the  unipolar  or  bipolar  reference  electrogram  recorded  from  the 
right  ventricular  epicardium  and  referred  to  the  onset  of  QRS.  A  character¬ 
istic  point  on  a  reference  electrogram  is  defined  as  a  distinct  peak,  notch, 
or  other  easily  identified  mark  that  remains  unchanged  over  the  entire  dura¬ 
tion  required  to  collect  electrograms  everywhere  on  the  epicardium.  The  time 
assigned  to  the  characteristic  point  on  the  reference  electrogram  is  measured 
from  the  onset  of  QRS,  usually  determined  from  the  onset  of  the  unipolar 
epicardial  lead  or  from  the  earliest  activation  recorded  from  intramural 
needle  electrodes  inserted  into  the  septum.  The  time  at  which  the  peak 
potential  occurs  for  electrograms  1  through  12  (Figure  A-3)  is  measured  with 
respect  to  a  characteristic  point  on  the  reference  lead  and  is  assigned  a 
negative  sign  if  it  occurs  before  the  characteristic  point  and  a  positive 
sign  if  it  occurs  after.  This  time  increment  is  added  algebraically  to  the 
characteristic  time  yielding  values  of  tn  measured  from  QRS  onset. 


Mathematical  Analysis 

The  principal  advantage  gained  by  arranging  the  electrogram  recording 
wires  into  a  geometric  array  is  that  a  mathematical  analysis  of  the  electro¬ 
gram  timing  will  yield  indirect  information  on  the  wavefront  properties. 
Information  about  the  epicardial  breakthrough  time  of  ventricular  excitation 
and  the  direction  of  propagation  over  the  epicardium  can  be  extracted  from 
computations  performed  on  12  electrograms.  Local  knowledge  of  the  time  and 
direction  of  the  excitatory  process  on  the  epicardium  is  calculated  for  each 
placement  of  the  clock  array.  Roving  the  clock  array  over  the  entire  epicar¬ 
dial  surface,  maintaining  a  constant  reference  electrogram,  produces  a  global 
picture  of  the  excitatory  process  from  which  the  ventricular  excitation  can 
be  inferred. 

The  relationship  between  the  electrogram  time  and  position  around  the 
clock  can  be  expressed  as  a  parametric  equation  over  the  surface  of  the  clock 
array.  A  mathematical  procedure  is  invoked  to  determine  the  best  estimate  of 
the  parameters.  Estimated  values  of  the  parameters  are  then  used  in  the 
original  equation  to  calculate  the  residuals  to  determine  the  validity  of  a 
single  wavefront  model  having  uniform  velocity.  When  the  residual  error  is 
small  compared  to  the  experimental  sources  of  error,  the  parametric  approxi¬ 
mation  provides  a  convincing  estimate  of  the  epicardial  breakthrough  time  and 
direction  of  movement. 

Epicardial  mapping  is  essentially  constructed  from  the  position  of  the 
clock  electrode  and  the  times  tn  contained  in  Figure  A-3.  On  some  portions 
of  the  epicardial  surface,  such  as  where  two  fronts  are  approaching  from 
different  directions,  the  best  way  to  construct  the  activation  sequence  is  to 
draw  isochronous  lines  over  the  clock  surface.  However,  for  the  major  por¬ 
tions  of  the  epicardial  surface  where  there  is  a  single  activation  front,  one 
would  like  to  analyze  the  activation  sequence  in  a  manner  that  minimizes  the 
influence  of  variability  in  the  local  times  tn  at  each  electrode.  The  varia¬ 
bility  in  the  time  tn  makes  it  difficult  to  construct  a  smooth  epicardial 
activation  sequence  by  interpolating  for  isochronous  locations  and  connecting 
them  with  straight  lines.  The  object  in  epicardial  mapping  is  to  obtain  a 


macroscopic  estimate  for  modeling  the  activation  over  a  finite  region  of  the 
epicardium  rather  than  activity  which  is  localized  to  a  point  on  the  surface. 
For  this  reason,  we  chose  to  analyze  the  times,  tn,  in  a  manner  that  utilizes 
an  algorithm  in  which  the  parameters  are  estimated  in  a  least  squares  sense. 
It  is  the  average  time  and  direction  over  the  surface  of  the  clock  electrode 
that  is  being  sought  to  construct  the  activation  sequence. 

If  the  ventricular  activation  front  propagates  with  a  uniform  conduction 
velocity  under  the  surface  of  the  clock  electrode  array,  then  the  times  tn 
will  bear  a  sinusoidal  relationship  across  the  face  of  the  clock  electrode. 
A  general  equation  for  the  time  t„  is  given  in  the  following  parametric  form: 

tn  =  T  sin  (w„  +  <t> )  +  B  (A-2) 

where  T  is  the  peak  time  measured  from  the  baseline  time  B.  wn  is  the  angle 
or  the  electrode  pair  measured  from  the  12  o'clock  position.  wn  =  30°x  n, 
where  n  =  1,2,,,, ,,12.  <t>  is  the  azimuth  or  phase  angle  which  defines  the 
direction  of  motion  across  the  clock  surface.  This  angle  <t>  i s  referenced  to 
the  12  o'clock  position.  B  is  the  baseline  or  mean  time  at  the  center  of  the 
clock  electrode.  Once  the  times  tn  are  recorded  from  the  clock  electrode, 
the  parameters  T,  <t>  ,  and  B  are  estimated  by  non-linear  least  squares.  Resi¬ 
duals  are  obtained  from  initial  estimates  of  the  parameters  given  by  the 
equation: 

ri~  *n—  TjSin(a)n  +  <l>j)  —  Bj  (A-3) 

Improved  values  of  the  parameters  T,  <J>  ,  and  B  are  obtained  from  the  data 
pairs  (tn,«n),  and  the  initial  estimates  T|  , Oj ,  B| ,  by  minimizing  the  resi¬ 
duals  r,  .  The  improvements  are  accomplished  by  a  differential  correction 
technique  based  on  least  squares  to  minimize  the  values  of  rj .  This  tech¬ 
nique  requires  the  initial  estimates  of  the  parameters  to  be  sufficiently 
close  to  the  actual  values  to  insure  rapid  convergence.  The  differential- 
correction  technique  obtains  updated  estimates  from  a  1 inear-Taylor  series 
expansion  of  function  (A-2)  around  the  parameters  T,<1>,  and  B. 


As  an  example  of  the  application,  consider  the  electrogram  recordings 
contained  in  Figure  A-3.  These  recordings  were  taken  from  the  left  ventricu¬ 
lar  wall  near  the  apex  of  a  normal  baboon  heart.  Each  epicardial  electrogram 
in  Figure  A-3  corresponds  to  a  clock  pair  of  recording  electrodes.  The 
number  associated  with  each  clock  pair  is  listed  in  column  one  of  Table  A-l. 
Column  2  contains  the  corresponding  time  for  which  the  electrogram  attains  a 
maximum  voltage.  This  time  is  measured  from  QRS  onset  determined  from  the 
unipolar  reference  electrogram  identified  in  Figure  A-3  as  channel  14.  Arri¬ 
val  times  at  the  clock  array  range  from  the  earliest  of  19  ms  for  pairs  7  and 
8.  This  lapsed  time  of  12  ms  across  the  1  cm  diameter  of  the  array  charac¬ 
terizes  the  phase  relationship  for  the  ventricular  activation  front  intersec¬ 
ting  the  epicardial  surface  in  the  direction  between  1:30  and  7:30  of  the 
clock  orientation.  These  raw  time  points,  listed  in  column  2,  appear  as 
circles  marking  the  vertical  axis  of  the  graph  in  Figure  A-4. 


TABLE  Arl.  SOLUTION  OF  EQUATION  (A-2)  FOR  THE  ESTIMATED  ACTIVATION  TIMES  FROM 
A  SET  OF  ELECTROGRAM  TIMING  tft . 
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Figure  A-4.  Output  data. 

(Above)  Activation  front  across  the  epicardial  surface  under  the  clock 
electrode. 

(Below)  Lower  curve  shows  visually  the  agreement  between  the  reconstructed 
curve  and  the  original  data. 
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Initial  estimates  T.  ,  <t.  ,  and  B;  are  calculated  from  the  times  t  given 
in  Table  A-l.  ...  " 


;the  average  of  all  12  values 


n 

Tj  =  1/2  (tmax  -  tmjn)  ;l/2  the  time  between  the  earliest  and  latest 
times  on  the  clock  surface. 


*i  =  30°  x  nmax  ;nmax  is  the  number  designating  the  electrode  pair 
which  has  the  latest  time. 


Figure  A-4  contains  two  forms  of  the  output  data.  The  upper  drawing 
indicates  the  activation  front  across  the  epicardial  surface  under  the  clock 
electrode.  The  direction  is  given  by  the  convergent  value.  The  2  ms  iso¬ 
chronous  lines  are  calculated  as  Tf  sin(oi  +  $f  )  +  Bf,  where  the  subscript  f 
refers  to  the  final  convergent  values.  The  lower  curve  shows  visually  the 
agreement  between  the  reconstructed  curve  and  the  original  data. 

Following  initialization  with  the  above  values,  the  parameter  estimation 
procedure  produced  convergent  values  of: 

Bf  =  25.2  msec 

Tf  =  6.1  msec 

4>f  =  222  degrees 

Inserting  these  convergent  values  of  B,  ,  T, ,  and  *f  into  equation  (A-2) 
produced  the  fitted  or  smoothed  values  of  tp  listed  in  column  3  of  Table  A-l. 
The  residuals  for  the  fitting  process  are  listed  in  column  4  for  each 
electrode  arrival  time. 


A  standard  deviation  of  Sv  =  .7  msec  is  obtained  from  the  list  of  resi¬ 
duals.  This  value  of  Sv  is  small  compared  to  the  time  required  for  transit 
across  the  electrode  array,  i.e.,  12.1  msec  for  transit  via  the  fitted 
values. 


A  relative  error  of  6%  (.7/12.1)  is  considered  small,  when  the  proposed 
use  for  the  data  is  in  modeling  efforts.  Comparing  Sv  to  the  transit  time  is 
a  reasonable  measure  of  whether  the  activation  is  a  single  wavefront  of 
uniform  velocity. 

From  the  analysis  of  these  12  electrograms  it  can  be  concluded  that  the 
epicardial  surface  at  the  point  where  this  array  was  positioned,  is  a  single 
wavefront  centered  about  Bf  =  25.2  msec  from  QRS  onset  and  is  propagating 
with  a  phase  angle  of  <j>(  =  222  degrees  relative  to  12  o'clock  and  with  nearly 
uniform  velocity. 


Animal  Model  and  Oata  Collection 


Intramural  activation  in  baboons  was  selected  as  a  basis  for  examining 
the  use  of  the  intramural  and  clock  electrodes.  These  animals  were  normal 


adult  males  of  the  papio  anubus  subspecies  and  weighed  25/30  kg.  Animals 
were  anesthetized  at  3/4%  to  1%  halothane  during  the  open-chested  recording 
period  which  lasted  approximately  30  minutes. 

Both  the  intramural  needle  and  clock  electrode  array  data  were  recorded 
for  later  comparison.  During  both  recording  sessions  the  right  ventricular 
reference  electrodes  were  recorded  simultaneously  with  the  corresponding 
clock  or  intramural  needle.  These  references  were  obtained  from  unipolar  and 
bipolar  electrodes  with  a  1-mm  separation  secured  to  the  right  ventricular 
surface  near  the  outflow  tract  of  the  right  ventricle.  The  time  of  the  peaks 
on  the  electrograms  of  both  the  clock  and  intramural  waveforms  was  measured 
from  the  same  peak  on  one  of  the  reference  electrodes.  Two  reference  elec¬ 
trodes  were  used  so  there  would  be  a  back-up  in  case  one  of  the  reference 
electrodes  was  to  become  dislodged  or  altered  during  the  course  of  the  entire 
data  collection  process.  It  also  provides  protection  against  changes  that 
might  be  orthogonal  to  a  single  electrode.  The  clock  electrode  electrograms 
were  collected  first  from  52  epicardial  sites.  These  sites  were  pre-selected 
and  located  on  a  p i c tor i al  view  of  the  epicardial  surface.  Approximately  10 
cardiac  cycles  were  recorded  from  each  site.  The  collection  process  advanced 
sequentially  across  the  52  epicardial  sites  in  a  continuous  recording  mode. 
Lapsed  time  for  the  entire  epicardial  data  collection  process  was  6  minutes. 
As  soon  as  the  epicardial  data  collection  was  completed,  the  intramural 
collection  started.  The  same  52  epicardial  sites  were  sequentially  mapped 
with  the  12-point  intramyocardial  needle  electrodes.  An  intramural  needle 
electrode  was  impaled  into  the  ventricle  at  the  same  location  the  epicardial 
clock  had  assumed  previously.  Several  moments  were  allowed  for  the  data  to 
stabilize  after  the  insertion  of  the  intramural  electrode  before  electrograms 
were  recorded  on  magnetic  tape.  Electrograms  from  the  septum  were  recorded 
from  5  sites  with  a  20-cm  long  intramural  needle  following  the  completion  of 
the  52  epicardial  sites. 


Intramural  sequence 

A  single  cross-section  of  the  intramural  activation  was  chosen  to  demon¬ 
strate  the  clock  array  and  compare  it  with  the  conventional  intramural  needle 
construction  of  the  activation  sequence.  The  complete  activation  sequence 
would  be  constructed  by  repeating  the  analysis  for  each  cross-section  of  the 
heart.  This  particular  cross  section  (Figure  A-5)  is  perpendicular  to  a  line 
from  the  apex  to  aorta  and  is  at  a  level  3  cm  above  the  apex.  The  activation 
sequence  depicted  in  Figure  A-5  was  constructed  from  12  ventricular  and  2 
septal  insertion  sites.  Timing  was  with  respect  to  the  onset  of  QRS  which 
was  measured  from  the  reference  electrogram.  Maximum  amplitude  of  each 
intramural  electrogram  marked  its  position  along  the  electrode  and  therefore 
its  time  relative  to  the  reference  electrogram. 


Phase  Velocity 

Extrapolating  an  intramural  activation  sequence  from  the  clock  array 
requires  the  approximation  of  a  uniform  conduction  velocity  for  the  excita¬ 
tory  process  within  the  myocardium.  In  the  case  of  the  baboon,  there  is  no 
estimated  value  available  for  both  open  and  closed  chest.  An  estimated  value 
was  obtained  by  impaling  the  20-cm  long,  intramural  needle  electrode 
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percutaneously  into  the  free  wall  of  the  left  ventricle  just  prior  to  opening 
the  chest.  The  electrograms  were  recorded  from  the  intramural  needle  simul¬ 
taneously  with  two  conventional  limb  leads.  As  soon  as  the  chest  was  opened, 
the  entrance  wound  from  the  intramural  needle  was  located  and  the  needle 
electrode  was  re-inserted.  Electrograms  were  recorded  again,  and  the  timing 
information  extracted  as  shown  in  Figure  A-6.  Conformity  between  open  and 
closed  chested  values  of  the  slope  of  this  line  (phase  velocity)  is  within 
the  uncertainty  of  the  measurements  and  hence  no  distinction  can  be  made 
between  the  apparent  differences  in  the  data.  Averaging  the  time  data  yields 
a  phase  velocity  of  0.32  mm/ms.  At  this  location  in  the  free  wall  of  the 
left  ventricle,  the  angle  of  incidence  of  the  propagating  front  was  estimated 
from  the  clock  array  to  be  8  degrees.  This  small  angle  of  incidence  indi¬ 
cates  a  phase  velocity  nearly  equal  to  the  conduction  velocity  at  this  ven¬ 
tricular  location.  In  subsequent  calculations  involving  the  clock  array  a 
value  of  0.32  mm/ms  is  taken  as  the  conduction  velocity  within  the 
myocardium. 


Extrapolated  Sequence 

Electrograms  recorded  with  the  clock  array  from  the  epicardial  sites  were 
analyzed  numerically  by  the  procedure  just  outlined.  The  12  epicardial  sites 
lying  on  the  same  cross-section  of  the  heart  illustrated  in  Figure  A-5  for 
the  intramural  activation  were  used  to  extrapolate  a  sequence  of  activation 
from  the  clock  electrode.  This  extrapolated  sequence  appears  in  Figure  A- 7 
where  the  value  of  the  uniform  conduction  velocity  was  taken  as  0.32  mm/ms. 
Table  A-2  lists  the  information  calculated  from  the  clock  array  which  is 
needed  to  estimate  the  ventricular  sequence,  mean  time  at  the  center  of  the 
clock  array,  and  the  azimuth  and  polar  angles  depicted  in  Figure  A-8.  The 
intramural  times,  although  not  used  in  the  extrapolation  process,  are  listed 
for  comparison  with  the  mean  epicardial  clock  times. 

Projecting  the  conduction  velocity  onto  the  plane  of  the  cross-section 
determines  the  phase  velocity  in  this  plane.  When  the  cross-section  is  taken 
normal  to  the  epicardial  surface,  the  phase  velocity  calculation  becomes 
Vph  =  u/cos  *  . 

Using  a  time  increment  of  5  ms  between  isochronous  lines,  the  displace¬ 
ment  increments  are  calculated  as  dx  =  5  u/cos  for  each  electrode  site. 
Figure  A-7  depicts  the  activation  at  5-ms  increments  and  is  to  be  compared 
with  the  intramural  activation  of  Figure  A-5. 

Comparing  of  the  mean  clock  times  and  the  corresponding  outer  intramural 
needle  times  indicates  a  discrepancy  in  the  two  epicardial  maps.  The  time 
uncertainty  between  these  epicardial  lines  is  not  considered  large  (2-4  ms) 
and  is  related  to  the  uncertainty  in  the  location  on  the  epicardium  of  the 
needle  and  clock  array.  When  the  polar  angle  is  large,  the  uncertainty 
produced  by  a  small  displacement  in  the  electrode  position  is  large;  i.e., 
where  the  propagation  velocity  across  the  epicardial  surface  is  slow,  the 
positioning  of  the  electrode  becomes  crucial.  The  purpose  for  comparing 
these  two  sequences  is  to  show  the  sequence  determined  from  two  independent 
methods  of  construction,  one  of  which  makes  use  of  a  modeling  simplification. 


TABLE  A-2.  ACTIVATION  SEQUENCE  EXTRAPOLATED  FROM  12  EP I  CARD  I AL  SITES  USING 
THE  CLOCK  ELECTRODE  ARRAY 
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namely  uniform  conduction  velocity.  On  a  macroscopic  level  of  observation, 
events  occurring  in  the  myocardial  activation  sequence  could  be  expected  to 
manifest  themselves  in  body  surface  ECGs.  Modeling  the  ECGs  of  a  normal 
baboon  with  a  propagation  program  utilizing  a  uniform  conduction  velocity 
would  not  produce  a  significant  difference  in  ECGs  generated  by  considering 
microscopic  variations  in  conduction  velocity.  The  extrapolated  ventricular 
sequence  of  Figure  A-7  is  sufficiently  similar  to  the  intramural  needle 
method  for  a  depth  of  10  mm  to  be  a  useful  approximation  of  the  wavefront  for 
a  forward  simulation  of  ECGs. 


WAVE  PROPAGATION  MODEL 


Heart  Geometry 

Since  those  electrical  events  that  are  simulated  by  the  forward  model  are 
confined  to  the  cardiac  muscle,  some  consideration  must  be  given  to  the  shape 
of  the  myocardial  mass.  The  human  heart,  specifically  the  right  and  left 
ventricles,  is  asymmetric  in  shape  and  therefore  models  such  as  spheres  and 
ellipsoids  are  unsatisfactory  representations  for  a  rigorous  simulation  of 
cardiac  electrical  events.  A  necessary  first  step  in  putting  together  a 
forward  model  was  to  obtain  an  accurate  human  ventricular  geometry,  represen¬ 
tative  of  normal  adult  males.  Such  a  ventricular  geometry  was  obtained  from 
a  42-year-old,  white,  male,  auto-accident  victim. 

Heart  fixing  and  sectioning  was  performed  by  first  filling  the  left,  and 
then  the  right,  ventricular  cavities  with  warm  formalin  agar-agar  at  the 
respective  normal  mean  diastolic  volumes.  The  heart  was  suspended  from  the 
root  in  a  square  container  of  warm  formalin  agar-agar  (see  Figure  A-9)  with 
the  A-V  groove  as  near  to  the  transverse  as  possible.  The  entire  mass  was 
allowed  to  cool  and  fix  over  night.  When  cooled  to  near  4°C,  the  block  of 
solid  gelatin  containing  the  heart  was  placed  in  a  rotary  blade  macrotome  for 
cutting  the  heart  breadloaf  fashion  at  2-mm  intervals  from  apex  to  base  (see 
Figure  A-10).  Total  heart  cross-sections  were  photographed. 

Drawings  with  the  estimated  normal  wall  boundaries  are  made  from  the 
photographs  of  each  2-mm  cross-section.  A  gross  anatomical  shape  of  the 
intact  myocardium  was  constructed  for  use  in  simulating  the  cel  1 -to-cel 1 
activation  pathway.  An  on-line  digitizing  table  was  used  to  transform  the 
anatomical  shape  into  a  suitable  digital  format.  After  digitization,  the  2- 
mm  cross-sections  parallel  to  the  A-V  groove  are  interpolated  to  1  mm  and 
this  becomes  the  anatomical  shape  and  resolution  of  the  ventricles 
incorporated  into  the  forward  model. 


Propagation  Model 

The  propagation  model  is  a  computer  program  that  simulates  the  myocardial 
activation  sequence.  This  simulation  includes  the  anatomical  and  physiologi¬ 
cal  factors  directly.  Parameters  such  as  intramural  conduction  velocity  and 
rate  of  Purkinje  activation  are  entered  as  first-order  approximat ions.  Reso¬ 
lution  and  accuracy  of  the  simulation  is  sufficient  to  produce  a  representa¬ 
tive  sequence  of  activation  for  a  particular  set  of  input  conditions,  such  as 
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heart  geometry  and  conduction  mechanism.  When  employed  as  a  mathematical 
model  of  the  heart  generator,  the  simulation  yields  data  on  the  area,  direc¬ 
tion,  and  locations  of  points  on  the  activation  front  at  each  instant  of 
time. 

Our  approach  to  the  problem  of  simulating  the  pathway  of  ventricular 
depolarization  was  to  develop  a  numerical  method  analogous  to  Huygen's  geome¬ 
trical  method  of  wavefront  construction  in  a  three-dimensional  region.  Antic¬ 
ipating  a  vast  amount  of  data  and  calculations,  the  numerical  method  had  to 
conform  to  the  restrictions  of  present  computer  configurations,  namely, 
memory  capacity  and  cycle  time.  Briefly,  Huygen's  method  consists  of  approx¬ 
imating  the  propagating  wavefront  by  a  large  number  of  spherical  wavelets 
(see  Figure  A-ll) . 

The  area  of  the  wavefront  is  taken  to  be  the  tangent  planes  connecting 
the  spherical  wavelets.  Wavelets  for  each  successive  time  increment  are 
determined  by  the  product  of  the  velocity  of  propagation  and  time  increment. 
This  distance  is  the  radius  of  the  spherical  wavelets  whose  centers  lie  on 
the  previous  wavefront.  To  simulate  a  method  such  as  this  with  a  resolution 
that  is  physiologically  meaningful  requires  a  computer  with  a  large  memory 
capacity,  since  the  location,  area,  and  direction  of  some  portion  of  each 
wavelet  must  be  computed  and  preserved  as  future  output  data.  The  digital 
computer  program  which  we  are  currently  running  is  essentially  a  replication 
of  Huygen's  method  with  isotropic  and  homogeneous  propagation  velocity.  The 
only  limitation  on  the  resolution  of  the  simulation  is  imposed  by  the  core 
capacity  of  the  computer.  This  comes  about  by  the  method  used  to  define  the 
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Figure  A-ll.  Spherical  wavelets. 
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regions  in  which  the  wavefront  propagates.  We  define  these  regions  by  their 
external  boundary  measurements  on  a  three-dimensional  Cartesian  coordinate 
system.  The  present  65K  RAM  machine  places  an  upper  limit  to  the  resolution 
of  1  mm  on  the  coordinate  system.  With  a  1-mm  coordinate  system  the  wavelets 
have  a  radius  of  1/2  mm  and  the  velocity-time  increment  is  1  mm.  We  confine 
the  ventricular  mass  to  lie  within  a  cube  which  is  100  mm  on  a  side.  To 
resolve  the  pathway  of  ventricular  excitation  to  within  1  mm  requires  the 
storage  of  1  million  coordinate  locations  in  memory. 

All  of  the  intersections,  of  coordinate  axes,  contained  within  the  exter¬ 
nal  boundaries  of  an  active  region  are  represented  by  "on"  bits  in  memory. 
Regions  where  no  activation  front  is  to  be  generated  are  represented  by  "off" 
bits.  As  the  activation  front  propagates  radially  outward  through  the  active 
region,  the  "on"  bits  are  switched  to  the  "off"  condition  until  the 
completion  of  all  activation. 

Wavefront  generation  is  initialized  by  specifying  the  coordinates  of  an 
arbitrary  set  of  points.  At  any  time  during  the  generation  of  the  wavefront, 
any  number  of  new  points  can  be  introduced  as  additional  sites  from  which 
wavefront  is  generated.  The  logic  employed  in  generating  the  propagating 
front  across  a  Cartesian  coordinate  system  is  based  on  finding  the  "on"  bits 
enclosed  on  the  surface  of  a  sphere  of  radius  1,  2,  or  3  mm.  After  each 
three  time  increments  the  process  is  started  over  using  the  last  "on"  bits  as 
the  center  for  the  next  three  spheres.  This  scheme  produces  a  maximum  propa¬ 
gation  time  error  along  the  diagonals  to  the  principal  planes  of  (2^2)/3  or 
less  than  5%. 

As  the  front  moves  out  across  the  intersections  of  the  coordinate  system 
in  a  radial  manner,  the  wavelets  are  constructed  at  each  intersection  for 
which  the  corresponding  bit  is  in  the  "on"  configuration.  Rather  than  com¬ 
pute  each  tangent  plane  connecting  each  spherical  wavelet,  some  portion  of 
the  surface  area  of  each  wavelet  is  used  to  approximate  the  tangent  plane. 
This  approximation  is  taken  to  be  the  projection  of  the  wavelet's  surface 
area  onto  one  of  the  principal  equatorial  planes.  The  resolution  of  the 
wavefront  is  the  minimum  increment  of  area  contributed  by  a  single  wavefront. 
This  minimum  increment  is  the  surface  area  of  one  octant  of  a  sphere  having  a 
radius  of  1/2  mm  that  is  "/8  mm.  From  the  unit  vectors  lying  in  the  coor¬ 
dinate  planes  and  located  at  the  center  of  the  wavelet,  the  resultant  direc¬ 
tion  of  the  wavelet  is  computed  by  taking  the  vector  sums  of  the  unit  vectors 
which  terminate  on  the  boundary  of  the  wavelet.  At  the  present  time  we  have 
chosen  as  output  the  total  area  of  the  wavefront  and  the  resultant  direction 
of  all  wavelets.  This  output  data  is  presented  in  the  form  of  vector  projec¬ 
tions  onto  the  principal  planes  of  the  coordinate  system.  Locations  of  the 
wavelets  at  each  instant  of  time  are  preserved  so  that  the  pathway  of  the 
activation  front  can  be  displayed. 


Computer  Program 

The  program  requires,  including  the  system  support  routines,  approximate¬ 
ly  65,000  locations  of  memory.  The  execution  time  varies  depending  on  the 
initial  condition.  Both  Fortran  IV  and  assembly  language  are  used.  A  major 
part  of  the  required  memory  is  used  to  store  the  geometry  of  the  object 
through  which  the  wave  will  propagate.  This  array  requires  130,000  bytes,  and 


is  dimensioned  100,  100,  13.  Using  each  bit  to  represent  a  location,  we  end 
up  with  a  cube  100  locations  on  a  side.  All  references  to  the  geometry  are 
as  if  there  were  100,  100,  100  words  in  the  array.  Assuming  the  indexing  to 
be  in  the  order  X,  Y,  Z,  the  location  (bit)  is  reached  by  using  the  value  of 
X  and  Y,  and  dividing  Z  by  8,  if  the  quotient  is  zero,  then  the  absolute 
location  is  the  8th  bit  of  that  byte.  If  the  quotient  is  not  zero,  add  one 
to  get  the  byte  and  the  remainder  is  the  bit  in  that  byte,  or  absolute 
location. 

Initially,  the  data  specifying  the  geometry  is  read  in  the  form  of  cards 
or  tape.  All  of  the  bits  in  the  cube  not  specified  by  this  geometry  are 
turned  off  or  set  to  zero.  This  establishes  the  criteria  that  the  front  may 
move  only  to  those  locations  containing  a  one  bit.  When  the  bit  is  zero, 
movement  to  that  location  is  prevented.  Next  comes  the  input  to  modify  the 
geometry.  Again  this  can  be  either  on  cards  or  tape.  This  input  allows  the 
reference  geometry  to  remain  fixed,  while  modifying  it  in  the  computer  from 
run  to  run. 

Starting  points  are  specified  immediately  following  the  geometry  section. 
These  points  specify  the  times  and  locations  for  starting  and  turning  on 
points  other  than  those  turned  on  as  the  front  moves  through  the  geometry. 
At  time  zero,  all  points  associated  with  that  time  are  taken  as  the  initial 
wavefront. 

The  first  step  in  the  simulation  is  to  compute,  considering  each  point  on 
the  front  individually,  each  incremental  area  and  direction.  The  six  bounding 
points  are  looked  at  for  each  point  on  the  front.  If  fewer  than  two  are  on, 
then  there  can  be  no  area.  If  five  or  more  are  on,  the  incremental  area  and 
direction  are  computed  and  stored.  If  there  are  more  than  two  points,  but 
fewer  than  five,  a  check  is  made  to  insure  that  the  points  do  not  lie  in  a 
plane,  allowing  the  front  to  degenerate  to  a  line.  Based  on  the  results  of 
this  test,  the  total  area  bounded  by  the  front  and  its  direction  is  computed 
from  the  incremental  values  and  the  results  stored. 

The  next  step  is  to  zero  out  the  points  on  the  old  front  and  compute  the 
coordinates  of  those  points  to  represent  the  new  front.  The  procedure  here 
is  to  assume  that  each  initial  point  is  the  center  of  a  sphere.  After  the 
first  iteration,  T  =  1,  the  sphere  will  have  a  radius  of  1.  All  of  the 
points  that  lie  within  a  radius  of  1  and  that  are  "on"  will  be  included  as  a 
point  on  the  new  front.  The  points  selected  are  stored  and  zeroed  out  and  a 
new  point  and  sphere  are  considered.  This  procedure  continues  until  all  of 
the  spheres  have  been  scanned.  All  of  the  new  points  are  returned  to  their 
"on"  status.  The  same  procedure  is  followed  after  the  completion  of  the 
second  iteration.  After  the  third  iteration,  the  same  procedure  for  selec¬ 
tion  of  the  points  on  the  new  front  is  performed,  except  these  new  points  are 
stored  as  initial  points  on  new  spheres  of  radius  one  for  the  next  iteration. 
A  radius  of  three  is  the  maximum  a  sphere  can  have.  This  is  also  necessary, 
as  a  radius  of  three  is  a  minimum  to  insure  that  no  points  are  missed  as  the 
front  moves  and  so  that  the  front  will  move  in  a  uniform  manner.  This,  in 
essence,  says  that  after  every  three  iterations,  the  simulation  is  restarted 
with  a  new  set  of  initial  values. 

The  final  step  in  the  completion  of  each  iteration  is  to  check  to  see  if 
any  additional  points  are  to  be  turned  on.  If  there  are,  then  the  parameters 


of  these  points  are  read  from  the  input  source  and  added  to  points  considered 
part  of  the  current  front.  The  simulation  program  continues  to  iterate  in 
the  above  manner  until  the  number  of  points  for  the  new  front  goes  to  zero. 


Program  Testing 

Accuracy  of  the  simulation  technique  was  studied  by  comparison  of  the 
results  obtained  from  propagation  within  an  isotropic  and  homogeneous  sphere. 
The  analytical  expression  for  the  wavefront  within  the  sphere  can  be  derived 
readily  when  the  initial  point  of  activation  is  on  the  surface  of  the  sphere. 
For  a  sphere  of  diameter  50  mm  and  the  initial  activation  point  on  the 
surface,  surface  area  of  the  wavefront  of  radius  r  enclosed  by  the  sphere  is 
given  by  the  equation: 


A  _  2tt  r2  [i  -  r/50] 


( A-4) 


Figure  A-12  shows  a  plot  of  equation  (A-4)  superimposed  on  an  equatorial 
plane  of  the  sphere.  A  smooth,  unbroken  curve  is  obtained  from  equation  (A- 
4)  and  represents  the  area  of  the  front  of  radius  r  corresponding  to  the 
smooth  outer  circle  which  defines  the  exterior  of  the  sphere.  The  jagged 
lines  bounding  the  sphere  from  the  outside  represent  the  approximation  of  the 
sphere  to  within  1  mm.  This  approximation  is  taken  as  the  geometry  for  the 
simulation  which  also  has  its  initial  activation  point  on  the  exterior  sur¬ 
face.  As  shown  in  Figure  A-12,  the  broken  line  indicates  the  surface  area  of 
the  wavefront  calculated  by  the  simulation  program  at  integer  values  of  the 
radius  r. 

Comparing  the  two  curves  during  the  early  portion  of  the  wavefront  devel¬ 
opment,  it  is  clear  that  the  area  computed  from  the  simulation  program  is 
larger  than  the  area  computed  from  equation  (A-4)  for  increasing  values  of  r, 
partly  because  the  approximation  to  the  sphere  is  larger  than  the  sphere 
itself.  The  error  in  this  approximation  is  zero  at  the  equatorial  plane 
which  is  perpendicular  to  the  direction  of  propagation.  From  equation  (A-4) 
the  value  of  the  radius  r  for  the  wavefront  as  it  passes  through  this  equa¬ 
torial  plane  is  given  as  33-1/3  mm.  This  occurs  at  the  extremum  for  the  area 
A  given  by  equation  (A-4).  The  true  area  of  the  wavefront  enclosed  by  the 
sphere  at  r  =  33-1/2  mm  is  2326  mm2,  whereas  the  simulation  program  computes 
a  value  of  2238  mm2,  or  approximately  4%  error. 

Around  a  radius  of  45  mm  the  disagreement  is  largest  between  equation  (A- 
4)  and  the  simulated  wavefront  area.  The  major  portion  of  this  error  is  due 
to  the  propagation  time  lag  of  5%  along  the  diagonals  to  the  principal 
planes.  This  time  lag  causes  the  wavefront  to  be  larger  than  it  should  be 
whenever  the  amount  of  active  region  is  decreasing  with  increasing  distance 
from  the  initial  activation  point.  In  general,  this  error  is  a  property  of 
the  geometry  and  is  indeterministic.  In  addition,  some  portion  is 
contributed  by  the  oversized  approximation  to  the  geometry  of  the  sphere. 

Since  the  diameter  of  the  sphere  is  50  mm  and  the  time-velocity  increment 
of  the  simulation  is  taken  as  1  mm,  there  were  50  time  points  at  which  the 
wavefront  area,  position,  and  direction  were  specified.  No  time  lag  was 
present  in  the  total  propagation  interval  because  the  terminal  portions  to  be 
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Figure  A- 12.  Plot  of  equation  (A-4)  superimposed  on  an  equatorial  plane  of 
the  sphere. 

activated  did  not  lie  along  the  diagonals  of  the  Cartesian  coordinate  planes. 
Direction  cosine  numbers  were  computed  for  the  resultant  wavefront  at  each  of 
the  50  points.  Since  the  initial  activation  point  was  chosen  to  lie  along 
one  of  the  coordinate  axes  (viz.,  the  Y  axis  ),  it  was  anticipated  that  the 
direction  cosine  number  would  be  (0,1,0)  respectively  for  the  (X,Y,Z)  axes.. 
To  within  three  significant  figures,  the  direction  cosine  numbers  were 
(0.000,  0.999,  0.000)  for  all  50  time  points. 

Another  analytical  case  considered  was  the  rectangular  parallelepiped. 
This  case  was  chosen  so  as  to  determine  how  well  the  simulated  wavefront 
would  approximate  a  plane  wave.  If  the  activation  is  initialized  at  the 
center  of  one  face  of  the  parallelepiped  and  allowed  to  propagate  down,  the 
front  is  very  nearly  a  plane  wave.  The  dimensions  in  this  case  were 
11  mm  x  50  mm  for  the  rectangular  paral lelepiped.  The  activation  simulation 
program  required  exactly  50  steps  to  reach  the  opposite  end  of  the  parallele¬ 
piped  and  the  direction  cosine  numbers  always  remained  (0.000,  0.000,  1.000). 
Near  the  end  of  the  simulation  the  area  of  the  wavefront  was  113.04  mm  which 
wnen  compared  to  the  actual  121.00  mm, results  in  an  error  of  approxima te I y 


Simulation  of  the  Purkinje  network  is  carried  out  separately  from  the 
activation  in  the  ventricular  muscle.  The  computer  program  is  identical  in 
Doth  cases,  but  the  initial  points  representing  the  conduction  bundles  and 
the  points  on  the  endocardial  surface  representing  the  Purkinje  network  are 


simulated  as  a  separate  phase  with  a  conduction  velocity  three  times  larger 
than  the  cell-to-cell  conduction  velocity  in  the  ventricular  wall.  Results 
from  the  first  phase  of  simulation  are  used  as  input  data  to  drive  the  second 
program  which  generates  the  activation  sequence  throughout  the  entire  cardiac 
geometry. 


Normal  Activation 


The  simulation  begins  by  selecting  the  input  parameters  to  phase  one  of 
the  computer  program.  This  first  phase  simulates  the  Purkinje  propagation 
over  the  endocardial  surface.  Only  the  coordinate  locations  on  the  endocar¬ 
dium  representing  Purkinje  are  contained  in  this  portion  of  the  program.  The 
effects  of  the  conduction  bundles  are  introduced  by  initiating  Purkinje 
geometry  at  three  endocardial  areas  in  the  left  ventricle  synchronously, 
followed  by  two  endocardial  areas  of  the  right  ventricle  and  right  septal 
surfac°.  Specifically,  the  input  parameters  are  five  arrays  of  coordinate 
points  on  the  endocardial  surface  from  which  the  propagating  wavefront  ori¬ 
ginates.  To  simulate  a  normal  ventricular  activation  sequence  in  the  human, 
we  chose  input  parameters  based  on  the  experimental  observations  of  Durrer  et 
al.  (A-4).  With  multi-point,  bipolar  electrodes,  they  mapped  the  time  course 
and  instantaneous  distribution  of  the  excitatory  process  on  seven  isolated 
human  hearts  having  no  history  of  cardiac  disease.  We  initiated  the  Purkinje 
network  of  this  simulation,  using  their  observations  that  endocardial  excita¬ 
tion  starts  at  five  sites,  three  of  which  are  in  the  left  ventricle:  high 
superior  paraseptal  area,  septal  anterior,  low  inferior  paraseptal  area,  and 
soon  followed  by  the  inner  wall  of  the  right  ventricle  and  the  septal  surface 
of  the  right  ventricle.  We  assumed  the  absence  of  any  Purkinje  penetration 
into  the  myocardium  and  set  the  Purkinje  conduction  velocity  at  120  cm  per 
second.  The  Purkinje  network  comprised  all  the  points  on  the  endocardial 
surface  except  for  a  region  of  the  base  of  the  left  septal  surface. 

The  digital  computer  program  simulation  of  the  activation  of  the  Purkinje 
network  during  phase  one  above,  produces  a  time  sequence  of  points  across  the 
endocardial  surface.  This  time  sequence  of  points  is  used  to  initiate  the 
simulation  of  the  ventricular  depolarization  in  phase  two  at  the  rate  of 
120  cm  per  second.  Simultaneously  the  myocardial  activation  is  propagating 
at  a  velocity  of  40  cm  a  second  isotropically,  and  homogeneously.  This 
produces  activation  fronts  which  are  essentially  from  endocardium  to  epicar- 
dium.  The  computer  program  then  continues  to  execute  the  activation  sequence 
according  to  Huygen's  method  until  all  the  points  of  the  cardiac  geometry 
have  been  traversed. 

During  each  iteration  of  the  computer  program,  the  wavefront  propagates  a 
distance  of  1  mm,  the  grid  size  upon  which  the  cardiac  geometry  is  resolved. 
At  the  conclusion  of  each  iteration,  the  program  calculates  the  total  wave- 
front  area  and  the  vectorial  area,  including  the  vector  location  of  the 
numerous  spherical  wavelets  with  a  radius  0.5  mm,  which  comprises  the  excita¬ 
tion  wavefront.  Selecting  the  ventricular  wavefront  velocity  of  40  cm  per 
second  yields  2.5  ms  of  time  between  iterations  and  therefore  isochronous 
wavefronts  separated  by  2.5  ms.  A  discussion  of  this  particular  effort  to 
simulate  normal  human  ventricular  activation  can  be  found  in  reference  ( A- 5 ) . 
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Two  views  of  the  activation  sequence  are  presented  in  Figures  A-13  and  A- 
14,  a  horizontal  and  a  vertical  perspective.  Isochronous  time  lines  at  5  msec 
intervals  represent  the  location  of  the  activation  front.  Conduction  velo¬ 
city  of  40  cm/sec  yields  a  time  increment  of  2.5  msec  for  each  1-mm  iteration 
step  in  the  activation  program.  Therefore,  the  isochronous  lines  correspond 
to  every  other  iteration  of  output  data  from  the  program. 

Only  the  gross  initializing  effects  of  the  right  and  left  conductiun 
bundles  are  reported  in  these  results.  Figures  A-13  and  A- 1 4  show  the  gross 
effects  of  the  bundles.  The  mid-left  septal  surface  and  right  septal  surface 
near  the  apex  are  activated  over  areas  approximating  the  regions  shown  in  the 
activation  sequence  map  in  humans  reported  by  Durrer.  Also,  the  effects  of 
the  3  to  1  ratio  between  Purkinje  and  ventricular  wall  velocities  can  be  seen 
as  a  rapid  movement  along  the  endocardial  surface  from  apex  to  base  and  from 
septum  to  left  free  wall.  Subendocardial  Purkinje  penetration  was  not 
included  in  this  preliminary  simulation. 

Qualitatively  at  least,  the  simulated  activation  sequence  map  corresponds 
to  observed  activation  in  normal  humans.  General  features  such  as  the  septal 
cancellation,  movement  from  apex  to  base,  from  inside  outward  are  visible  in 
Figures  A-13  and  A-14. 


EQUIVALENT  CARDIAC  GENERATOR 


Discussion 

Associated  with  the  excitatory  pathway  in  the  ventricles  is  the  myocar¬ 
dial  cell  membrane  action  potential.  The  currents  responsible  for  the  ECG  on 
the  torso  surface  arise  from  the  EMFs  generated  by  the  cellular  action  poten¬ 
tials.  Figure  A- 1 5  is  an  example  of  a  myocardial  action  potential,  obtained 
from  Burgess.  Figure  A- 1 6  is  a  plot  of  a  mathematical  model  of  the  action 
potential  (A-P)  from  a  continuous  function  of  time.  Specific  details  of  the 
mathematical  model  will  be  discussed  later  under  the  subject  of  equivalent 
generators.  Equivalent  cardiac  generators  depend  on  the  shape  of  A-Ps  and 
are  important  to  the  simulation  studies  conducted  with  the  forward  model. 
Varying  the  shape  of  the  A-Ps  within  regions  of  the  ventricles  will  result  in 
body  surface  ECGs  via  the  simulation  that  reflects  these  changes. 

A  convention  exists  in  the  literature  for  discussing  the  various  phases 
of  an  action  potential.  As  will  become  clear  in  the  discussion  of  the  A-P 
equation,  the  model  parameters  are  directly  related  to  the  phases  of  the  A-P, 
hence  the  definitions  used  in  the  forward  model  will  be  stated  at  this  time. 
In  reference  to  Figure  A-16,  the  upstroke  of  the  cardiac  A-P  is  labeled  phase 
0  and  is  the  sudden  depolarization  of  the  cell  following  a  stimulus  that 
raises  the  cell-membrane  potential  to  threshold,  by  way  of  the  gap  junctions 
of  the  intercalated  discs  separating  the  interiors  of  adjacent  cells.  Phase 
1  is  a  brief  phase  of  rapid  repolarization  that  occurs  within  a  few  milli¬ 
seconds  following  depolarization.  The  prolonged  duration  of  phase  2  is 
characterized  by  a  plateau  of  nearly  stabilized  membrane  potential.  Phase  3 
marks  the  final  repolarization  process  of  the  cell  membrane.  Following 
repolarization,  the  cell  remains  in  its  phase  4  resting  state.  Equivalent 


generators  for  the  recovery  current  are  obtained  from  equations  of  the  A-P, 
the  parameters  in  these  equations  determine  the  shapes  of  the  various  phases. 

The  stimulus  that  invokes  a  large  increase  in  membrane  permeability  to 
sodium  (A-6,A-7)  in  one  region  of  the  heart  initiates  activity  that  spreads 
in  all  directions  at  a  velocity  depending  on  the  cable  and  excitable  proper¬ 
ties  of  the  cells.  Activity  spreads  through  cardiac  muscle  from  cell-to-cell 
and  the  tissue  acts  electrically  as  a  single  unit,  hence  the  term  "functional 
syncytium"  is  used  to  characterize  the  pathway  of  depolarization.  In  syncy¬ 
tial  tissue,  depolarization  spreads  out  from  the  originating  point  as  a 
wavefront  in  three-dimensional  space.  The  propagation  program,  one  of  the 
programs  comprising  the  forward  model,  is  a  simulation  of  this  functional 
syncytium  and  the  three-dimensional  wavefront  of  depol arization. 

An  equivalent  cardiac  source  (ECS),  for  the  purpose  of  the  forward  model, 
is  a  mathematical  expression  for  the  cardiac  EMF  which  replaces  the  signifi¬ 
cant  physical  properties  of  the  actual  cardiac  source.  Since  the  forward 
model  is  essentially  a  numerical  process,  the  description  of  the  source 
impressing  current  into  the  torso  volume  conductor  must  eventually  be  reduced 
to  a  numerical  process  likewise. 

Selecting  an  ECS  is  guided  by  three  principal  considerations:  (1)  EQU I  VA¬ 
LENCE — the  mathematical-physical  equation  serving  to  replicate  the  current 
source  must  be  equivalent  in  the  sense  that  the  equation  of  continuity  of 
current  density  is  not  violated;  (2)  REPRESENTATION--the  time  varying  spatial 
distribution  of  the  actual  current  source  must  be  approximated  to  some  speci¬ 
fied  order  that  establishes  how  well  the  equivalent  source  represents  the 
actual  current  source;  (3)  MATHEMATICAL  F0RM--the  mathematical  equation  for 
the  ECS  should  lend  itself  to  physical  interpretation,  i.e.,  the  coefficients 
and  variables  should  relate  to  physical  phenomena  and  be  quantitatively 
determined  by  experimentation.  Many  trade-offs  and  approximations  are  re¬ 
quired  that  compromise  the  above  objectives  and  result  in  ECSs  that  must  be 
evaluated  with  respect  to  a  set  of  objectives  for  the  forward  simulation. 

One  of  the  simplest  ECS  models  that  has  been  considered  for  a  suitable 
description  of  the  electrical  source  in  the  heart  is  the  single  dipole.  The 
single  dipole  is  the  algebraic  sum  of  all  the  elemental  bipolar  current 
sources  in  the  heart  at  any  one  instant.  A  fixed  position  dipole  with  three 
time-varying  orthogonal  components  is  the  source  model  employed  in  vectorcar¬ 
diography.  A  major  problem  with  the  fixed  position  dipole  model  is  the  lack 
of  information  on  the  spatial  distribution  of  the  elemental  bipolar  sources. 
Spatial  dimensions  within  the  heart  are  not  small  when  compared  to  the  dis¬ 
tances  to  the  chest  surface.  Although  the  single  fixed  dipole  model  is 
equivalent  and  has  a  tractable  mathematical  form,  the  representation  of  the 
actual  source  produces  a  poor  approximation.  A  single  moving  dipole  model  is 
less  objectionable  but  still  not  a  sufficiently  good  approximation. 

Equivalence  and  representation  of  the  ECS  is  well  defined  for  the  multi¬ 
pole  expansion  model.  Since  the  multipole  model  is  attained  from  a  series 
expansion  of  the  Coulomb  potential  equation  in  a  Taylor  series,  equivalence 
is  possible  by  eliminating  the  monopole  term  and  the  ECS  representation  by 
the  series  is  complete  in  a  mathematical  sense.  Only  the  number  of  terms  in 
the  expansion  determines  the  order  of  approximation,  which  can  be  determined 
analytically,  and  compared  to  the  measurement  system  errors.  The  difficulty 
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with  a  multipole  expansion  model  is  the  mathematical  form  is  undesirable 
since  there  is  no  anatomical  or  electrophysiological  correlates  to  the  coef¬ 
ficients  in  the  Taylor  series  expansion,  which  are  the  multipole  moments.  A 
model  such  as  this  replaces  empirical  body  surface  ECG  methods  with  empirical 
heart  multipole  methods. 

More  recently,  forward  and  inverse  model  studies  have  been  performed  with 
the  multiple  dipole  model  as  the  equivalent  cardiac  source.  As  the  name 
implies,  multiple  (finite  number)  dipoles  are  distributed  strategical ly  at 
significant  sites  within  the  myocardium.  In  this  manner  of  approximation, 
each  dipole  represents  a  region  of  myocardium  that  has  small  dimensions  when 
compared  to  the  distance  to  the  chest  surface.  A  multiple  dipole  model  obeys 
the  equivalence  criteria  and  has  a  mathematical  form  that  is  tractable  from  a 
numerical  analysis  perspective.  Representation  by  the  multiple  dipole  source 
depends  on  numerous  heart-torso  variables,  but  is  essentially  determined  by 
the  number  of  dipoles  and  their  respective  locations.  There  are  no  signifi¬ 
cant  constraints  on  the  number  of  dipoles  that  can  be  programmed  into  the 
forward  model. 


Bipolar  Model 


Since  the  propagation  model  has  a  spatial  resolution  of  1  mm,  the  mathe¬ 
matical  model  of  choice  for  the  ECS _wi 1 1  be  a  bipolar  source  of  EMF  with  a  1- 
mm  separation.  Given  two  points  r,  and  fj  in  the  myocardium  separated  by 
1  mm,  the  mathematical  form  of  the  bipolar  source  is  derived  from  the  poten¬ 
tial  at  these  two  points.  Consider  the  potential  at  point  r,  as  Vj  and 
s  imi  larly  for  r ■  .  V; 


A  bipolar  source  vector  b  is  defined  for  a  finite  separation  of  d,  taken  as 
1  mm  in  the  propagation  model.  Vector  b  has  the  following  form: 

V  [vn>-vn»]  7 


where  the  direction  of  bu 


is  from  r-  to  r(  and  is  given  by 

—  r  •  _  T 

d  = 


_ Ll 

!ri  -  ril 


An  equivalent  dipole  source  for  a  segment  of  the  myocardium  is  constructed 
from  the  vector  sum  of  the  bipolar  sources  within  the  nth  segment,  i.e., 


P(i  ( rn 


ij 


V. 
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where  the  point  r,  is  summed  over  six  adjacent  points,  repeated  for  all  M 
points  in  the  segment.  In  the  forward  model  simulation,  the  number  of  points 
M  in  each  segment  varies  with  the  size  and  shape  of  each  segment,  which  in 


turn  depends  on  the  arbitrary  subdivision  of  the  myocardium.  With  more 
segments  in  the  heart  and  fewer  points  in  each  segment,  an  increase  in 
resolution  is  achievable  (hence  improved  representation  of  cardiac  electrical 
events)  at  the  expense  of  computer  resources  and  effort  to  specify  the  ini¬ 
tial  conditions  to  the  simulation  program.  Each  composite  dipole  moment  Pn 
of  the  nlh  segment  serves  as  both  depolarization  and  repolarization  equiva¬ 
lent  cardiac  source,  since  no  distinction  has  been  made  between  depolariza¬ 
tion  and  repolarization  potentials.  The  vector  position  rn  of  the  equivalent 
dipole  is  the  mean  value  of  the  m  triplets  (xi,yi,zi)  of  coordinate  points 
within  the  nfh  segment.  It  is  this  value,  rn,fo r  which  the  boundary  value 
problem  of  a  unit  source  computes  the  transfer  impedances  to  the  torso 
surface  electrode  sites. 

Depolarization  and  repolarization  of  the  cardiac  cells  produces  a  system 
of  charges  and  currents.  A  simply  connected  surface  can  be  specified  that 
completely  surrounds  the  cardiac  cells.  The  fixed  volume  Q  bounded  by  the 
surface  r  contains  all  the  charges  and  currents  comprising  the  cardiac  source 
as  a  function  of  time.  It  can  be  shown  that  the  distributed  cardiac  system 
of  charges  and  currents  can  be  represented  by  a  time-varying  dipole.  In 
addition,  the  volume  Q  can  be  subdivided  into  a  set  of  smaller  volume  ele¬ 
ments  and  a  time-varying  dipole  can  be  defined  for  the  charges  and  currents 
inside  each  element.  Since  the  equivalent  cardiac  sources  in  the  forward 
model  are  taken  as  dipole  moments,  it  is  essential  that  the  cardiac  system  of 
charges  and  currents  be  rigorously  represented  as  time-varying  vector 
quantities. 

The  dipole  moment  of  the  charge-current  distribution  is  defined  by 


p  -  f  r  P  dn 


(A-5) 


where  r  is  the  position  vector  from  a  fixed  origin  and  p  is  the  volume  charge 
density  inside  the  fixed  volumes.  Consider  J  to  be  the  current  density 
inside  Q  .  Then  as  an  application  of  Green's  theorem  on  the  fields  j  and  r  , 
the  following  identity  results: 


/  J  dfi  r  (f  r  ( J.U'!  dr  -  f  r  OIV  j  dfi 
A  Jv 


( A- 6 ) 


j  vanishes  on  the  surface  r  because  it  is  contained  completely  inside  n  . 
Equation  (A-6)  reduces  to 


The  right-hand  side  of  equation  (A-7)  can  be  rewritten,  using  the  equation  of 
continuity  for  charges 


at 


+  DIV  J  -  O 


and  substituting for  Div  J. 
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(A-8) 


The  partial  derivative  of  p  with  respect  to  time  becomes  the  ordinary  time 
derivative  of  the  integral  whe  brought  outside  the  integral  sign. 


fn  J  dn  =  4—  f  r  pdn 

•'fi  dt  Jo 


(A-9 ) 


By  definition,  the  right-hand  integral  is  the  definition  of  a  dipole  moment 
of  a  charge-current  distribution  given  by  equation  (A-5). 


_  dP 


_  J  dn  = 
n  dt 


(A-10) 


This  derivation  establishes  the  validity  of  using  time-varying  dipoles  as 
the  equivalent  cardiac  source  for  the  distributed  charge-current  within  the 


heart.  Equation  (A-10)  is  the  basis  for  representing  the  electrical  activity 
in  the  heart  by  an  equivalent  set  of  time-varying  bipolar  sources.  Briefly 
stated,  the  depolarization  and  repolarization  charge-current  distribution  J 
is  generated  by  the  cell  membrane  EMF  which  is  confined  to  the  myocardial 
volumes,  i.e.,  the  space  occupied  by  ventricular  myocardium.  Regardless  of 
the  spatial  distribution  or  time  variation  of  the  cellular  EMFs  producing  the 
charge-current  J  ,  the  resulting  source  (the  intejral  of  J  over  n  )  jemains 
equivalent  to  a  time-varying  dipole-vector  sourcedP/dt  •  The  vector p  ,  when 
expressed  as  a  combination  of  linearly  independent  vectors,  represents  the 
multiple-dipole,  equivalent  cardiac  source. 

The  forward  simulation  of  ECGs  uses  an  equivalent  cardiac  source  that  is 
a  depolarization-repolarization  generator  defined  as  a  multiple  bipolar 
source  (no  limiting  process  invoked  as  with  the  dipole)  of  EMF  within  the 
myocardium  which  is  a  function  of  time  measured  from  the  onset  of  depolariza¬ 
tion.  Equivalent  current  is  expressed  by  an  equation  that  is  a  continuous 
function  of  time  and  the  shape  of  the  equation  is  determined  by  the  values  of 
six  parameters.  A  set  of  six  parameters  are  specified  for  each  coordinate 
intersection  in  three  dimensions  within  the  myocardium.  Given  that  the  heart 


geometry  is  specified  on  a  1-mm  grid  system,  th 
extends  over  a  distance  of  1  mm  between  sourc 
intersection  within  the  myocardial  geometry  is 
current.  Figure  A- 1 7  depicts  a  typical  interse 
dial  geometry.  An  equivalent  bipolar  moment  is 
EMF  at  the  intersecting  point  constructed  fro 
intersection  points  surrounding  some  arbitrary 
the  reference  point  or  location  for  the  equival 
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Similarly  for  Py  and  Pz 


py  =CV2-  V0]  d  +  [V,-  Vo]  d 
pz  -  [  - V0]  d  +  [V.-V0]d 


Hence  the  vector  moment  of  the  equivalent  bipolar  source  formed  from  Px  ,  Py , 
and  P z  is 

P  =  Px'T'  +  Pyt  +  P2t 


where  i,  j,  and  k  are  unit  vectors  along  the  x,  y,  and  z  axes,  respectively. 
The  above  expression  contains  the  spatial  property  for  the  bipolar  source. 
As  a  function  of  time,  P(i,j,k,t)  varies  implicitly  as  the  voltages  Vo(t), 
Vi (t) , , , ,V6(t)  vary,  while  the  separation  d  is  held  constant.  A  mathematical 
equation  expressing  the  voltage  at  each  coordinate  intersection  as  a  contin¬ 
uous  function  of  time  is  employed  to  compute  each  Vn ,  and  therefore 
P(i»j  >k,t) . 

Action  Potential  Profile 

Figure  A-18  is  a  graph  of  the  mathematical  equation  Vn  for  a  typical  set 
of  parameters. 

[V,  ;  t  <  T, 


V„  = 


h  +  a'(t-T1)  ?  T,  <  t  <  t2 

V2  ♦  cosf.i7(t +  io  +  Vt  ea3(t~^-a* 

l  +  eaj(t  -T2)  -  a3 


Time  T-|  is  obtained  from  the  propagation  simulation  program;  it  is  the 
time  relative  to  an  arbitrary  reference  (usually  the  onset  of  QRS)  at  which 
the  activation  front  crosses  the  point  (i,j,k). 

ai  =  (V2  -  V1  )/dt  dt  =  2.5  ms,  phase  0 

a2  =  refractory  rate,  phase  3 

a3  =  functional  refractory  time,  phase  3 

V-,  =  resting  potential  difference 

V2  =  depolarization  potential  difference 


500  ms 


Figure  A-18.  Action  potential  waveforms  as  generated  by  Vp  on  p.  114  .  In 
the  top  set  of  waveforms,  the  parameter  a3  (the  functional 
refractory  time,  phase  3)  was  changed.  In  the  middle  set  of 
waveforms  az,  the  slope  of  phase  3,  was  varied.  In  the  bottom 
panel  both  were  varied  as  shown. 
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n  of  the  myocardial  geometry  is  characterized  by  a  set 
An  equivalent  depolarization-repolarization  source 
entire  myocardial  mass  is  a  matter  of  practical  choice 
ts  to  resolve  the  spatial  distribution  of  EMF  and  as 
te  intersections,  of  which  there  are  approximately 
the  equivalent  cardiac  source  consists  of  20  segments: 
le  and  septum,  and  8  for  the  right  ventricle.  Each  of 
left  ventricle  and  septum  is  divided  into  zones.  The 
are  divided  into  3  zones  from  endocardium  to  epicar- 
segments  are  divided  likewise  into  4  zones,  and  the 
nto  5  zones.  In  summary,  the  segment-zone  subdivision 


TRANSFER  IMPEDANCES 


Boundary  Value  Problem 


The  total  bioelectric  current  density  J  within  the  volume  conductor 
contributed  by  two  components  given  by  the  following  vector  equation: 


J  =  J; 


( A- 1 1 ) 


J,  is  the  impressed  current  density  resulting  from  the  nonconservative  elec¬ 
tric  field  of  depolarization/repolarization  such  that  a  path  integral,  equa¬ 
tion  (A-12),  through  the  source  region  yields  the  imbedded  electromotive 
force. 


/Ei 


.dL  =  EMF 


(A-12) 


The  quantity  o  E  is  the  ohmic  current  density  present  in  the  volume  conductor 
with  conductivity  o  due  to  the  conservative  electric  field  E  produced  by  \t  . 

Applicability  of  Ohm's  law  to  the  electrocardiographic  forward  problem  is 
based  on  the  nominal  values  of  conductivity  and  permittivity  of  volume  con¬ 
ductor  tissue  and  the  observed  maximum  rate  of  change  of  electrical  events 
within  the  heart.  A  full  discussion  of  this  area  is  given  by  Plonsey  (A-8). 
He  points  out  that  nominal  values  of  conductivity  and  permittivity  of  tissue 
yield  a  sufficiently  short  time  constant  to  permit  the  forward  problem  to  be 
considered  quasi -static .  This  means  that  volume  charge  density  changes  or 
capacitive  effects  can  be  neglected.  Forward  solutions  calculated  indepen¬ 
dently  for  each  time  point  will  not  have  appreciable  time-delay  errors  caused 
by  the  volume  conductor. 

Setting  the  volume  charge  density  rate  of  change  9t  equal  to  zero 
results  in  a  stationary  (non-t ime-dependent )  partial  differential  equation. 
The  definition  of  current  and  the  law  of  conservation  of  charge  leads  to  the 
fol lowing  definition: 
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where  the  current  dq/dt  results  from  the  time  rate  of  change  of  volume 
current  density  contained  within  the  volume  F . 


Applying  Gauss'  theorem  to  equation  ( A- 1 3 )  results  in  the  equation  of 
conti nui ty: 


v.,  -  fr  =  0 


( A- 14 ) 


For  a  stationary  problem  3o/3t  =  0  and  equation  (A-14)  reduces  to  a  diver¬ 
gence  equation  for  the  total  current  density: 


V  •  J  =  0 


( A- 1 5 ) 


Instituting  for  the  total  current  density  from  equation  (A-ll),  a  nonhomo- 
:erecus  partial  differential  equation  appropriate  to  the  forward  problem  is 

;pr ■ ved. 


V.  E  =  -  v.  J: 


(A- 16) 


iauat'on  (A-16)  can  be  restated  in  terms  of  a  potential  since  the  conserva- 
::ve  * i e 1 d  E  is  derivable  from  the  gradient  of  a  scalar  potential,  say  V. 
’rere^ore,  equation  (A-16)  becomes 


V*  <7  v  y  =  v  •  J  j 


(A- 17) 


For  purposes  of  discussion  equation  ( A- 1 7 )  will  be  written  in  the  form  of  an 
operator  equation: 


=  f 


( A- 1 8 ) 


where  the  bioelectric  current  sources  impressed  by  the  generator,  7 • J ,  ,  are 
replaced  by  the  free  term  symbol  f  and  the  partial  differential  operation  on 
the  potential  function  V  is  replaced  by  the  operator  symbol^. 

Equation  ( A- 1 8 )  leads  to  a  boundary-value  problem  when  the  heart-volume 
conductor  is  imbedded  in  an  insulating  medium.  The  volume  conductor-to-a i r 
interface  imposes  a  boundary  condition  on  the  normal  component  of  the  current 


density.  For  the  torso  air  interface  where  o=  0,  continuity  of  the  normal 
component  of  current  density  requires  it  to  be  zero.  That  is 


J  (torso)  =1  (air)  =  0 

n  n 


where  the  subscript  n  refers  to  the  normal  to  the  surface  S,  and  from  Ohm's 
1  aw 


Jn  =  *  En  ~ 


“  —  <J 


av 

an 


This  yields  the  natural  boundary  condi 


=  0 

tion. 


an 


s 


=  0 


where  S  is  the  surface  bounding  the  volume  conductor  y  . 
boundary-val ue  problem  can  be  stated  as: 


The  resulting 


<E  V  =  f  in  * 


av 


a  n 


=  0 


(A- 19) 


It  is  well  known  from  operator  theory  in  partial  differential  equations  that 
a  solution  to  problem  (A- 19 )  exists  if  the  operatorj^is  positive  definite. 
This  existence  criteria  imposes  a  restriction  on  the  source  term  f,  namely, 


f 2  d'l'  < 


(A- 20) 


The  equality  of  sources  and  sinks  of  current  or  non-unipolar  nature  of 
the  heart  current  sources  insures  this  physically.  The  choices  of  a  parame¬ 
tric  mathematical  model  for  the  forward  problem  are  restricted  to  equations 
that  satisfy  conditions  ( A- 1 0 ) . 

Uniqueness  can  be  established  by  imposing  a  reference  condition: 


V  d*  = 


0 


( A- 21) 


Gelernter-Swihart  Method 


The  Gelernter  and  iwihart  method  is  a  numerical 
homogeneous  Toisson  equation  for  an  inhomogeneous 
divides  the  volume  conductor  into  regions  of  consta 


solution  to  a  piece-wise 
volume  conductor.  It 
nt  conductivity  and  then 
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calculates,  for  an  arbitrary  source  distribution,  the  potential  distribution 
on  an  irregularly  shaped  surface  bounded  by  an  insulator.  Each  homogeneous 
region  R  is  formulated  as  a  Neumann  problem  for  Poisson's  equation: 


v2V  =  f 


in  R 
on  C 


(A-22) 


where 

9  v  /  3  n 


R  is  a  bounded  connected  region  with  sufficiently  smooth  boundary  C  and 
denotes  the  normal  derivative  taken  in  the  outward  direction.  It 


follows  from  Green's  first  identity  that  f  and  g  cannot  be  chosen  indepen¬ 
dently,  but  must  satisfy  the  relation 


I'™ 


-  f  i  dC 
JC 


( A- 23) 


The  solution  to  equation  (A-23)  is  only  unique  up  to  an  additive  con¬ 
stant.  This  constant  is  usually  determined  by  a  normalization  process  over 
the  surface  such  that  the  integral 


/  V  dC  =  0 
yC 


(A-24) 


Boundary  conditions  are  of  two  types;  the  external  boundary  requires  the 
normal  component  of  the  current  density  to  be  zero,  whereas  the  internal 
interfaces  require  the  normal  component  to  be  refracted  in  proportion  to  the 
conductivity  on  each  side  of  the  interface. 


Gelernter-Swihart  Algorithm 

The  GS  procedure  for  computing  the  solution  to  the  Neumann  problem  of 
Poisson's  equation  is  quite  general  in  scope  and  is  applicable  to  the  forward 
ECG  problem  if  practical  limitations  are  imposed  on  the  number  of  internal 
regions  requiring  specification  of  geometry  and  conductivity,  although  no 
such  requirements  are  imposed  by  the  GS  method  itself.  Derivation  of  the 
method  considers  an  arbitrary  ensemble  of  generators  in  a  closed,  three- 
dimensional  region  containing  sub-regions  that  each  have  a  uniform  conduc¬ 
tivity.  The  region  is  required  to  be  embedded  in  an  insulating  medium. 
There  are  no  restrictions  other  than  closure  upon  the  configuration  of  the 
sub-regions. 

Solving  the  boundary  value  problem  follows  closely  the  actual  events 
which  bring  about  equilibrium  in  a  charge  distribution  for  a  physical  process 
having  the  same  properties  that  are  described  numerically  in  the  computer 
program. 

Consider  the  following  as  a  brief  discussion  of  the  Gelernter-Swihart 
formulation,  based  on  the  computer  program  rather  than  on  their  formal  treat¬ 
ment  (A-9).  A  dipole  of  known  location  and  size  is  contained  inside  the 
bounded  inhomogeneous  volume  conductor.  Initially,  the  charge  distribution 


induced  on  all  surface  area  elements  in  an  unbounded  medium  is  calculated. 
Gauss'  induction  law  is  employed: 


E(Xj )*a$j  cos (  *,)  =  q. 


where  E  is  the  electric  field  produced  by  the  dipole  at  the  surface  area 
element  ASj.  4> 4  is  the  angle  between  E  and  the  unit  normal  of  as*  .  qj  is 
the  induced  charge  on  the  area  element  ASj  .  Q0  is  a  functional 
representation  for  all  the  qj.  E  for  a  dipole  is  given  by: 


3(  p  •  x  )  x  p 


The  distribution  of  charge  Q0  on  the  boundaries  is  then  redistributed 
using  Coulomb's  law  for  the  electric  field  of  point  charges  and  Gauss'  induc¬ 
tion  until  the  final  charge  distribution  is  in  equilibrium.  This  is  accom¬ 
plished  numerically  by  combining  Coulomb's  and  Gauss'  laws  so  that  the  charge 
induced  on  area  element  As,  by  a  charge  at  Xj  can  be  expressed  by  the  solid 
angle  the  area  A$f  subtends  at  Xj ,  i.e.,  Qij, 


Aq(  -nu  qj 


For  the  inhomogeneous  regions,  n,j  is  multiplied  by  (o1  -o2)/(a1  +  a2), 
the  ratio  of  conductivities  at  the  interfaces.  This  algorithm  is  used  in  an 
iterative  manner  until  the  induced  charges  go  to  zero. 

Beginning  with  Q0  the  induced  charges  on  the  first  iteration  become 


a  Qi  =  ft  Qo 

a  Q2  =  QaQ., 

second  iteration: 


A  Qn  =fiAQn.1 

as  n  becomes  large  a  Qn->- 0.  When  convergence  has  been  attained,  the  potential 
distribution  is  calculated  from  the  final  charge  distribution  over  all  sub- 
regions  and  bounding  region. 

No  convergence  criteria  were  provided  for  the  GS  method  and  since  no 
exactly  calculated  result  was  available  for  comparison  purposes,  they  offered 
the  following  arbitrary  stopping  criteria  for  the  charge  iteration  algorithm: 


1.  The  difference  in  potential  between  the  maximum  and  minimum  points 
changes  by  less  than  an  acceptably  small  amount  from  its  value  at  the 
last  iteration. 

2.  The  normal  component  of  the  electric  field  at  each  surface  element 
has  reached  an  acceptably  small  percentage  of  the  maximum  normal  field 
produced  by  the  uncompensated  generator  (i.e.,  the  boundary  condition  at 
the  external  surface  is  acceptably  satisfied). 

3.  The  surface  potential  map  has  become  essentially  stable  (i.e.,  the 
equipotential  lines  and  maximum  and  minimum  points  have  not  changed  their 
position  noticeably  since  the  last  iteration). 

Acceptably  small  in  the  above  specifications  depends  upon  the  order  of 
approximation  required  of  the  solutions  in  a  particular  application. 

Rotes  of  convergence  were  estimated  from  solutions  of  simple  problems  to 
solutions  of  more  realistic  examples  and  ranged  from  4  to  40  iterations. 
High  conductivity  ratios  in  the  vicinity  of  the  generators  required  the  most 
iterations  before  the  stopping  criteria  were  met. 

Analytical  solutions  can  be  used  to  test  the  program  code  of  the  GS 
formulation;  however,  the  errors  attributed  to  boundary  and  sub-region  des- 
cretization,  stopping  criteria,  and  source  vector  position  for  the  arbitrary 
torso  configuration  are  in  doubt  without  a  known  solution  for  comparison. 
The  electrolytic  tank  with  immersed  bipolar  sources  provides  a  physical 
situation  for  which  a  GS  solution  can  be  computed  and  compared  to  a  measured 
potential  distribution.  The  torso  tank  used  to  test  the  computer  code  for 
mistakes  and  obtain  error  estimates  in  computed  solutions  is  a  life-sized 
Lucite  model  of  a  male  torso.  Bipolar  sources  are  constructed  and  calibrated 
with  a  constant  current  in  a  spherical  tank  before  positioning  them  in  the 
torso  tank.  A  potential  distribution  on  the  torso  tank  surface  is  measured 
from  64  uniformly  placed  electrodes  and  compared  with  a  corresponding  poten¬ 
tial  distribution  calculated  by  the  method  of  GS.  This  torso-shaped  electro¬ 
lytic  tank  test  bed  is  not  intended  to  be  a  physical  replication  of  the 
human;  its  only  purpose  is  to  provide  a  completely  and  accurately  specified 
problem  as  input  to  the  GS  algorithms,  for  which  a  measured  solution  exists. 


Transfer  Impedances  for  QRS  and  T  Wave  Simulation 

Simulation  of  the  T  wave  requires  that  the  motion  of  the  heart  during 
contraction  be  taken  into  consideration.  Shortly  after  the  ventricles  de¬ 
polarize  and  begin  repolarizing,  the  ventricles  begin  contracting  and  chang¬ 
ing  position  in  the  chest.  Since  the  equivalent  cardiac  generators  represen¬ 
ting  the  repolarization  EMF  must  remain  in  the  ventricular  myocardium,  their 
positions  must  follow  that  of  the  ventricular  wall.  Simulation  of  the  ECG 
during  the  T  wave  requires  transfer  impedances  from  the  repolarization  gen¬ 
erators  to  the  body  surface.  With  the  heart  moving,  the  repolarization 
generators  have  a  time-varying  position  during  contraction.  The  forward 
model  of  T  waves  takes  into  account  time-varying  transfer  impedances  by 
computing  3  sets  of  transfer  impedances  --one  at  diastole,  mid-systole,  and 
end-systole -- and  then  interpolates  between  these  3  sets  for  intermediate 
transfer  impedances  over  the  entire  contracting  period  of  the  heart. 


Data  on  the  location  and  motion  of  the  equivalent  cardiac  generators  is 
available  in  several  publications  (A-9-11),  where  Ingels  et  al.  measured  the 
midwall  myocardial  dynamics  in  24  intact  patients  by  radiographic  analysis  of 
surgically  implanted  markers.  They  implanted  tiny  coils  of  tantalum  wire 
into  the  left  ventricular  myocardium  in  24  patients  at  the  time  of  cardiac 
surgery.  Seven  coil  markers  were  placed  in  the  myocardium;  at  the  left 
ventricular  apex,  and  at  three  equidistant  points  from  apex  to  base  along 
both  the  anterolateral  and  inferior  margins  of  the  left  ventricle  in  a  pat¬ 
tern  outlining  the  ventricular  chamber  as  seen  in  a  30-degree  right  anterior 
oblique  projection.  Postoperatively  the  motions  of  the  myocardial  markers 
were  visualized  in  a  resting  single-plane  (30-degree  right  anterior  oblique) 
fluoroscopic  to  video  recording  at  30  frames  per  second.  The  positions  of 
all  myocardial  markers  for  each  frame  of  three  complete  cardiac  cycles  were 
used  to  calculate  the  positions  of  midwall  left  ventricular  segments. 

Figure  A-19  defines  the  length  and  radii  for  which  Ingels  et  al.  pub¬ 
lished  quantitative  data  over  the  entire  cardiac  cycle,  although  this  figure 
illustrates  the  definitions  transformed  into  the  ventricular  segments  of  the 
heart  model  utilized  by  the  simulation  and  shown  in  Figure  A-2.  Variations 
in  the  radii  defined  in  Figure  A-l  were  graphed  as  a  function  of  time  for  3 
cardiac  cycles.  The  values  of  the  radii  at  diastole,  mid,  and  end  systole 
are  taken  from  the  graphics  for  each  ventricular  segment  4  through  9  inclu¬ 
sive.  Radii  for  segments  10  to  12  and  1  to  3  are  taken  as  the  average  of  the 
two  segments  at  the  corresponding  positions  from  apex  to  base.  For  example, 
the  radius  of  segment  11  is  the  average  of  the  radii  for  5  and  8. 


Parameters  such  as  changing  radii,  longitudinal  axis  shortening,  and  the 
rotation  of  the  transverse  ventricular  diameters,  become  the  elements  of  a 
transformation  matrix  applied  to  the  equivalent  cardiac  source  locations  in 
the  simulated  heart  geometry.  After  performing  the  contraction  matrix  trans¬ 
formation  on  the  cardiac  source  locations,  the  results  are  transformed  again 
into  the  body  orientation  by  the  same  matrix  used  to  specify  all  heart 
geometry  information  in  the  chest.  These  transformations  yield  the  x,y,z 
coordinates  of  the  cardiac  sources  within  the  inhomogeneous  volume  conductor. 
Transfer  impedances  are  computed  as  the  solution  to  a  boundary  value  problem 
for  the  new  cardiac  source  locations.  Computing  these  solutions  at  diastole, 
mid  and  end  systole  produces  three  complete  sets  of  transfer  impedances.  All 
transfer  impedances  from  QRS  onset  to  the  end  of  the  T  wave  are  taken  from 
this  set  of  three  or  are  interpolated  from  them  at  the  appropriate  time 
increment. 


ECG  Lead  Systems 

ECG  tracings  and  vectorcardiograms  are  comprised  of  voltages  measured  as 
potential  differences  on  the  chest  surface.  Similarly,  the  forward  model 
produces  simulated  ECGs  and  vectors  that  are  computed  as  potential 
differences . 

Given  a  set  of  computed  potentials  over  the  chest  surface,  a  mathematical 
transformation  must  be  performed  on  these  potentials  to  produce  the  appro¬ 
priate  voltages.  Transforming  chest  potentials  into  ECG  and  vector  voltages 
requires  a  consideration  of  two  aspects  of  ECG  recording  equipment.  First 
the  voltages  are  determined  from  a  specific  assortment  of  potential 
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,  and  secondly  the  voltages  are  combined  by  weighting  factors 
determined  from  a  set  of  resistors  in  the  recording  equipment.  The  transfor¬ 
mation  of  chest  potentials  to  ECG  and  vector  leads  consists  of  a  matrix  array 
in  which  the  row  elements  are  the  weighting  factors  for  each  lead  and  the 
columns  correspond  to  the  various  chest  potential  sites. 


differences . 
determined  from 
ma 


Simulating  the  convention  12-lead  ECG,  McFee,  Frank,  Cube,  and  SVEC  III 
vectors  requires  38  potential  sites  in  the  chest  surface.  Potentials  from 
these  38  chest  sites  are  transformed  into  24  weighted  voltages  (12  leads  plus 


Lead  voltages  for  the  entire  ECG  are  represented  by  a  matrix  V(j,n), 
where  j  spans  the  rows  covering  the  chest  electrode  sites  and  n  specifies  the 
column  corresponding  to  each  instant  in  time.  The  matrix  V(j,n)  is  defined  as 
the  product  of  the  transformation  matrix  T(j,k)  and  the  potential 
distribution  matrix  S(k,n). 


V(j,n)  =  T(j,k)  *  S(k,n) 


Potential  distribution  S(k,n)is  a  matrix  array  computed  by  combining  the 
heart  current  sources  P(i,n)  and  the  unit  potential  distributions  U(k,n,i). 
Index  (i)  ranges  over  the  bipolar  sources  i  =  1  to  60.  Indices  (k)  and  (n) 
have  the  same  meaning  as  above,  k  =  1  to  38  electrode  sites,  n  =  1  to  N  ECG 
time  steps.  As  a  matrix  array,  S(k,n)  is  defined  by  the  following  equation: 


S(k,n)  =  ^Pu(k,n,i)  *  P(i,n) 
i 


Unit  distributions  U(k,n,i)  are  computed  as  solutions  to  the  boundary 
value  problem  for  unit  sources  imbedded  in  the  volume  conductor  at  locations 
specified  as  diastole,  mid-systole,  and  end-systole.  From  these  three  solu¬ 
tions,  a  set  of  interpolated  solutions  are  computed  for  each  time  point 
between  diastole  and  end-systole.  Let  any  time  point  during  the  ECG  be 
represented  by  the  subscript  n.  Hence,  the  ECG  lead  voltages  V(j,n)  become 


v(j,n)  =  T(j,10  *^u(k,n,i)  *  P(i,n) 

i 


To  summarize: 

1.  V(j,n)  is  the  matrix  array  consisting  of  the  conventional  12- 
lead  electrocardiogram,  and  the  vector  cardiograms  of  McFee,  Frank,  Cube,  and 
SVEC  III  positioned  in  the  rows  indexed  by  j  =  1  to  38.  Each  column  contains 
the  voltages  for  a  particular  time  point  of  which  there  are  a  total  of  N 
columns,  indexed  by  n  =  1  to  N  ECG  time  divisions. 

2.  T(j,k)  is  the  matrix  array  containing  the  electrode-lead 
weighting  factors  that  transform  a  chest  surface  potential  distribution  into 
ECG  lead  voltages. 

3.  P ( i , n )  is  the  matrix  array  whose  row  elements  are  the  bipolar 
source  moments  produced  by  the  simulation  of  heart  current  generators  (i  =  1 
to  60)  for  each  instant  of  time  (n),  the  column  index. 


4.  U(k,n,i)  is  the  matrix  array  having  row  elements  corresponding 
to  the  unit  potential  distribution  over  38  electrode  sites  for  unit  bipolar 
source  (i),  that  are  interpolated  from  diastole  to  end-systole  to  provide  the 
column  elements  corresponding  to  each  instant  of  time  (n).  The  array  has  a 
depth  (i)  covering  the  60  unit  sources. 

It  is  the  array  V(j,n)  that  is  displayed  on  output  as  the  ECGs  and  vectors 
for  depolarization  and  repolarization  of  the  ventricles. 
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